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1  Introduction 


Although  mammography  is  a  valuable  screening  tool  for  breast  cancer,  it  is  less  effective 
in  younger  women  (<  40  years),  usually  because  the  higher  density  of  their  breasts  can 
obscure  tumors.  While  the  incidence  of  cancer  in  younger  women  is  relatively  low,  it  is  still 
significant:  an  estimated  13,000  women  under  40  will  be  diagnosed  with  breast  cancer  this 
year.  Further,  cancers  in  this  age  group  tend  to  be  more  aggressive,  and  survival  rates  are 
lower.  In  light  of  these  observations,  there  is  need  for  an  effective  screening  technique  to 
complement  clinical  and  self  breast  exams.  In  addition,  the  lack  of  a  specific  early  detection 
screening  protocol  for  young  women  is  particularly  disquieting  to  those  at  elevated  risk. 

Elasticity  Imaging  (El)  [1]  is  a  technique  that  could  potentially  assume  this  role,  as  it 
relies  on  extracting  information  from  ultrasound  images  which  are  in  turn  unaffected  by  the 
denseness  of  the  breast.  The  application  of  El  to  breast  cancer  detection  (see  [2]  for  example) 
utilizes  the  fact  that  breast  tumors  tend  to  be  significantly  stiffer  than  the  surrounding 
tissue  [3],  and  can  be  easily  discerned  in  an  image  that  represents  the  spatial  variation  of 
stiffness  properties  [4],  Measurements  for  El  can  be  made  by  only  slightly  modifying  the 
protocol  for  a  typical  sonogram,  and  involve  recording  ultrasound  images  of  the  breast  as 
it  is  deformed,  and  registering  these  images  to  produce  a  displacement  field.  Using  this 
displacement  field,  and  an  assumed  linear  elastic  model  for  the  tissue,  an  inverse  problem 
is  solved  to  compute  the  spatial  variation  of  the  elastic  modulus.  Thus  El  yields  additional 
information  about  the  tissue  (elastic  modulus)  at  little  or  no  extra  clinical  cost. 

While  several  clinical  studies  are  currently  underway  to  assess  the  utility  of  El  as  a 
screening  tool,  we  believe  that  in  order  to  realize  its  true  diagnostic  potential  the  assumption 
of  linear  elasticity  within  El  must  be  replaced  by  a  more  realistic  non-linear  model.  In 
particular,  producing  images  of  material  parameters  E  (the  elastic  modulus  at  zero  strain) 
and  m  (the  degree  of  non-linearity  of  tissue)  for  a  non-linear  model  where  the  stress  (a) 
is  expressed  in  terms  of  the  strain  (e)  as,  cr  =  (. E/m)(eme  —  1),  will  allow  for  a  clearer 
differentiation  of  several  important  tissue  types.  For  example,  it  is  difficult  to  distinguish  a 
ductal  carcinoma  in-situ  from  a  phyllodes  tumor  based  on  E  alone  (since  the  ratio  is  about 
1.1),  while  it  should  be  relatively  easy  to  do  so  based  on  m  (the  ratio  is  about  2).  The 
opposite  holds  for  an  infiltrating  ductal  carcinoma  and  an  intraductal  papilloma,  where  the 
ratio  of  E  is  about  2.1,  and  the  ratio  of  m  is  0.93  (values  taken  from  [3]).  By  replacing  the 
linear  response  in  the  inverse  problem  in  El  with  a  more  informative  and  accurate  non-linear 
response  we  propose  to  generate  multiple ,  quantitative  images  which  provide  complimentary 
information  about  the  tissue,  instead  of  a  single  elastic  modulus  image. 

It  is  well  known  that  the  microvasculature  vascular  density  is  a  strong  predictor  of 
prognosis  in  cancerous  tissue.  It  can  be  show  that  this  quantity  is  directly  related  to  the 
microvascular  filtration  coefficient,  which  determines  the  rate  at  which  fluid  is  carried  away 
in  the  vascular  compartment  and  hence  the  rate  at  which  the  tissue  relaxes  to  an  applied 
load.  Recently,  a  biomechanical  model  which  accounts  for  this  fluid  transport  mechanism 
and  the  corresponding  relaxation  has  been  proposed  [5].  Using  this  model  we  now  wish  to 
verify  whether  the  strain  relaxation  of  tissue  due  to  fluid  transport  is  sufficiently  large  to 
be  measured  by  standard  clastography  techniques.  Our  hypothesis  is  that  if  this  spatio- 
temporal  relaxation  of  strain  can  indeed  be  measured  using  elastography,  then  an  inverse 
algorithm  can  be  designed  based  on  the  biomechanical  model  of  Netti  et  al.,  which  will  yield 
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the  distribution  of  the  microvascular  filtration  coefficient  in  the  tissue.  This  approach  may 
offer  another  exciting  avenue  for  the  application  of  elastography  to  cancer  detection. 


2  Body 

The  research  carried  out  under  this  grant  is  divided  into  two  components.  One  deals  with 
the  extension  of  elastography  to  non-linear  elasticity  imaging  and  the  other  aims  to  ascertain 
whether  elasticity  imaging  can  be  used  to  determine  the  microvascular  filtration  coefficient 
(MFC). 

2.1  Nonlinear  Elasticity  Imaging 

Onr  long  term  goal  is  to  develop  and  test  an  ultrasound-based  methodology  for  generating 
multiple-parameter,  high-contrast  elasticity  images  of  breast  tissue  for  improved  diagnosis 
and  detection  of  breast  cancer  in  young  women.  In  the  current  grant  we  have  obtained  proof 
of  concept  results  for  this  technique.  Or  results  are  described  in  detail  in  Appendix  A. 

2.2  MFC  feasibility  study 

Here  our  goal  was  to  determine  whether  changes  in  the  microvascular  filtration  coefficient 
in  cancerous  tissue  produce  detectible  changes  in  strain  relaxation  patterns,  which  may  be 
measured  using  standard  elastography  techniques.  This  research  paves  the  way  for  using 
elastography  to  determine  the  spatial  distribution  of  MFC.  The  results  of  this  research  are 
described  in  Appendix  B. 


3  Key  Research  Accomplishments 

1.  Non-linear  models  for  breast  tissue:  We  have  demonstrated  that  the  Veronda- 
Westmann  model  [6]  has  a  strain  energy  term  that  yields  the  experimentally  observed 
exponential  stress-strain  behavior  observed  in  most  breast  tissue  [3]. 

2.  Algorithms:  The  use  of  a  non-linear  model  adds  computational  complexity  to  the 
inverse  problem.  We  have  developed  and  implemented  a  quasi-Newton,  adjoint 
approach  to  solve  this  problem. 

3.  Test  problems:  Using  an  in-house  suite  of  computer  programs,  we  have  generated 
synthetic  displacement  data  (with  noise)  for  tissue  models  with  known  elastic 
parameters.  We  have  used  this  data  to  test  the  performance  of  the  inversion  algorithms, 
and  the  feasibility  of  the  overall  approach. 

4.  Poroelastic  model  for  tissue  deformation:  We  have  implemented  the  model  of  Netti 
et  al.  in  a  finite  element  program  in  order  to  simulate  the  relaxation  of  tissue  due  to 
vascular  drainage  and  interstitial  percolation. 
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5.  Feasibility  Study:  We  have  determined  that  vascular  drainage  (as  opposed  to  interstitial 
percolation)  is  the  main  mechanism  of  relaxation  in  most  vascular  tissue.  Through 
numerical  simulations  we  have  also  verified  that  the  increased  value  of  the  microvascular 
filtration  coefficient  in  cancerous  tissue  produces  changes  in  strain  relaxation  patterns, 
which  may  be  measured  by  standard  clastography  techniques. 


4  Reportable  Outcomes 

4.1  Manuscripts 

1.  Nonlinear  elasticity  imaging  I:  Formulation  and  computational  solution  for 
compressible  media,  by  Gokhale,  N.H.,  Oberai,  A. A.,  and  Barbone,  P.E.,  in 
preparation. 

2.  Nonlinear  elasticity  imaging  II:  Using  enhanced  strain  FEM  to  treat  nearly 
incompressible  media,  by  Gokhale,  N.H.,  Oberai,  A. A.,  and  Barbone,  P.E.,  in 
preparation. 

3.  Nonlinear  elasticity  imaging  III:  Inversions  in  nearly  incompressible  media  by  higher 
order  FEM,  by  Gokhale,  N.H.,  Oberai,  A. A.,  and  Barbone,  P.E.,  in  preparation. 

4.  Nonlinear  elasticity  imaging  in  tissue,  by  Gokhale,  N.H.,  Oberai,  A. A.,  and  Barbone, 
P.E.,  in  preparation. 

5.  Coupling  between  elastic  strain  and  interstitial  fluid  flow:  Ramifications  for  poroelastic 
imaging,  by  Leiderman,  R.,  Barbone,  P.E.,  Oberai,  A. A.,  and  Bamber,  J.  C.,  submitted 
to  Physics  in  Medicine  and  Biology. 

4.2  Presentations 

1.  Towards  the  early  detection  of  breast  cancer  in  young  women,  DOD  Era  of  Hope 
Conference,  Philadelphia,  June,  2005. 

2.  Biomechanical  Imaging:  Querying  Mechanical  Properties  of  Soft  Tissue,  In-vivo, 
Seminar  at  the  Department  of  Mechanical  Engineering,  Louisiana  State  University, 
Baton  Rouge,  Feb.  18,  2005. 

3.  Motion  and  Deformation  Measurement  from  Image  Sequences,  Advanced  Computation 
Seminar  Series,  Center  for  Computational  Science,  BU,  18  March  2005. 

4.  Inferring  Biomechanical  Properties  from  Quasistatic  Deformations:  An  introduction 
to  associated  Inverse  Problems.  Paul  Barbone,  Invited  Keynote  Tutorial,  Fourth 
International  Conference  on  the  Ultrasonic  Measurement  and  Imaging  of  Tissue 
Elasticity,  Austin  TX,  16  October  2005. 

5.  Nonlinear  Elasticity  Imaging,  Nachikct  Gokhale,  Assad  Oberai,  Paul  Barbone.  Eighth 
US  National  Congress  on  Computational  Mechanics,  Austin  TX.  25-27  July  2005. 
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6.  Progress  in  Biomechanical  Imaging,  Mike  Richards,  Nachiket  Gokhale,  Carlos  Rivas, 
Ricardo  Leiderman,  Assad  Oberai,  Paul  Barbone.  2nd  International  Conference  on 
Tumor  Progression  &  Therapeutic  Resistance.  18-20  September  2005. 

7.  Three  Dimensional  Ultrasound  Image  Registration  and  Shear  Elastic  Modulus 
Reconstruction,  Mike  Richards,  Nachiket  Gokhale,  Carlos  Rivas,  Ricardo  Leiderman, 
Assad  Oberai,  Paul  Barbone.  Fourth  International  Conference  on  the  LUtrasonic 
Measurement  and  Imaging  of  Tissue  Elasticity,  16-19  October  2005. 

8.  Progress  in  Biomechanical  Imaging.  Mike  Richards,  Nachiket  Gokhale,  Carlos  Rivas, 
Ricardo  Leiderman,  Assad  Oberai,  Paul  Barbone.  Fourth  International  Conference  on 
the  LUtrasonic  Measurement  and  Imaging  of  Tissue  Elasticity,  16-19  October  2005. 

9.  Progress  in  Quantitative  Biomechanical  Imaging.  Paul  Barbone,  Mike  Richards, 
Nachiket  Gokhale,  Carlos  Rivas,  Ricardo  Leiderman,  Assad  Oberai.  IMA  Workshop: 
Imaging  from  Wave  Propagation,  Institute  for  Mathematics  and  its  Application,  17-21 
October  2005. 

10.  Nonlinear  Elasticity  Imaging,  Nachiket  Gokhale,  Assad  Oberai,  Paul  Barbone,  ACES- 
IGERT  Review  Meeting,  29  September  2005. 

4.3  Degrees  Obtained 

Nachiket  Gokhale,  PhD.  in  Mechanical  Engineering  at  Boston  University,  2006.  Nachiket 

will  graduate  around  August  2006. 

4.4  Personnel  Supported 

1.  Dr.  Assad  Oberai  (PI) 

2.  Dr.  Ricardo  Leiderman  (Post-Doc) 

3.  Nachiket  Gokhale  (PhD.  Student) 

4.5  Funding  Applied  for  Based  on  Work  in  this  Award 

Title:  Biomechanical  Imaging,  PI:  P.E.  Barbone,  Boston  University,  Period: 

12/01/0611/30/11,  Source:  National  Institutes  of  Health  Total  Amount:  2,574,709. 

5  Conclusions 

Through  the  research  supported  by  this  award  we  have  developed,  implemented  and  tested 

(on  synthetic  data)  an  approach  to  determine  the  nonlinear  elastic  properties  of  breast  tissue. 

The  next  logical  step  will  be  to  test  this  approach  in  an  experimental  setting  using  tissue 

mimicking  phantoms  and  then  in  a  clinical  setting.  If  successful  in  a  clinical  setting  this 
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approach  will  substantially  increase  the  specificity  of  quasi-static  elastography  by  offering 
more  than  one  parameter  with  which  to  distinguish  different  tissue  types. 

In  addition  we  have  established  that  the  reported  variations  in  microvascular  filtration 
coefficient  (MFC),  which  is  a  strong  predictor  of  prognosis  in  cancerous  tissues,  lead  to 
significant  changes  in  strain  relaxation.  Further,  these  changes  in  strain  relaxation  are 
large  enough  so  that  they  may  be  detected  using  standard  ultrasound-based  elastography 
techniques.  This  motivates  the  development  of  an  imaging  modality  based  on  quantifying 
the  spatial  variation  of  MFC  which  could  play  a  significant  role  in  cancer  detection. 
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Appendix  A 


Iterative  solution  of  the  non-linear  inverse  elasticity 
problem 


Nachiket  H  Gokhale  §,  Paul  E  Barbone  and  Assad  A  Oberai  || 

Department  of  Aerospace  and  Mechanical  Engineering,  Boston  University,  Boston, 
MA  02215, USA 

Abstract.  We  discuss  and  solve  an  inverse  problem  in  non-linear  elasticity  imaging 
in  which  we  recover  spatial  distributions  of  hyperelastic  material  parameters  from 
measured  displacement  fields.  This  problem  has  applications  to  elasticity  imaging  of 
soft  tissue  because  the  strain  dependence  of  the  material  parameters  may  potentially 
be  used  to  differentiate  between  malignant  and  normal  tissues.  We  account  for  the 
geometric  and  material  non-linearity  of  the  tissues  by  assuming  a  known  hyperelastic 
model  for  the  soft  tissue.  We  formulate  the  problem  as  a  minimization  problem. 
The  cost  function  represents  the  difference  between  the  measured  and  predicted 
displacement  fields.  We  minimize  it  with  respect  to  the  spatial  distribution  of  material 
properties.  We  solve  the  minimization  problem  by  a  gradient  based  (quasi  Newton) 
optimization  approach.  We  calculate  the  gradient  efficiently  using  the  adjoint  method. 
We  present  numerical  examples  that  demonstrate  the  feasibility  and  efficiency  of  the 
approach,  and  compare  it  to  linear  elastic  reconstructions. 


1.  Introduction 

Elasticity  Imaging  (El)  is  an  emerging  medical  imaging  technique  in  which  images  of 
the  spatial  distribution  of  the  clastic  modulus  or  stiffness  of  soft  tissues  are  created. 
The  motivation  and  interest  for  persuing  this  goal  are  numerous.  For  example,  it  is  well 
known  that  pathologies  affect  the  mechanical  properties  of  soft  tissue.  This  is  evident 
in  breast,  prostate  and  other  tumors  presenting  as  hard  lumps,  in  fibrosis  which  is 
associated  with  a  diffuse  stiffening,  and  atherosclerosis,  which  means  literally,  hardening 
of  the  arteries.  There  is  also  interest  in  being  able  to  evaluate  the  normal  properties 
of  many  tissues  whose  function  is  primarily  mechanical,  such  as  lungs,  blood  vessels, 
muscles,  (including  cardiac  muscle),  and  cartilage.  Furthermore,  there  is  now  strong 
evidence  that  phenotypical  cell  behaviors  depend  relatively  strongly  on  the  mechanical 
properties  of  their  immediate  environments  [1].  Finally,  technological  applications  such 
as  automated  needle  insertion  [2]  or  surgical  planning  and  simulation  would  benefit  from 
patient  specific  a  priori  or  even  real  time  mechanical  characterization  of  soft  tissues. 
The  typical  procedure  followed  for  elasticity  imaging  is: 

§  To  whom  correspondence  should  be  addressed  (gokhalen@bu.edu) 

||  At  Rensselaer  Polytechnic  Institute,  Troy,  NY,  from  January  2006 


Iterative  solution  of  the  non-linear  inverse  elasticity  problem 
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(i)  Tissue  Deformation  and  image  acquisition:  The  soft-tissue  of  interest  is  imaged 
while  it  is  deforming.  The  source  of  the  deformation  in  general  may  be  external 
(e.g.  manual  palpation)  or  internal  (cardiac  motion).  The  deformation  may  be 
quasi-static  or  transient. 

(ii)  Image  Processing:  The  displacement  field  is  calculated  from  the  acquired  images, 
by  using  either  correlation  based  algorithms,  or  minimization  of  a  suitable  objective 
function  [3]. 

(iii)  Inverse  Problem  Solution:  The  spatial  distribution  of  the  material  properties  is 
calculated  from  the  known  (i.e.  measured)  displacement  field.  An  appropriate 
mathematical  model  for  the  tissue  deformation  is  selected  depending  upon  the 
spatio-temperal  scales  of  the  applied  deformation  and  the  imaging  system.  Given 
the  forward  mathematical  model  and  measured  deformations,  an  inverse  problem 
thus  follows.  In  this  paper,  an  efficient  formulation  based  on  the  adjoint  formulation 
is  used  to  solve  this  inverse  problem. 

Thus,  elasticity  imaging  involves  the  solution  of  the  following  inverse  problem:  Given 
the  displacement  field  in  an  elastic  body,  calculate  the  spatial  distribution  of  material 
properties  (stiffness).  Most  approaches  used  in  the  literature  to  image  the  stiffness  [4-10] 
assume  that  the  tissue  can  be  modeled  as  a  linear  incompressible  isotropic  clastic  solid. 
The  goal  then  is  to  recover  the  shear  elastic  modulus  distribution  of  the  soft  tissue. 
The  linear  elastic  model,  of  course,  can  accurately  predict  small  strain  deformations. 
There  are  two  expected  practical  advantages,  however,  to  using  large  strains  in  elasticity 
imaging.  One  is  that  the  signal-to-noise  ratio  of  the  measured  deformations  tends  to 
increase  with  deformation  magnitude.  In  this  case,  even  if  the  stress-strain  response  of 
the  tissue  remains  linear  over  the  range  of  strains  measured,  the  error  due  to  geometric 
nonlinearity  (by  this,  we  mean  the  neglect  of  the  nonlinear  terms  in  the  exact  definition 
of  strain  [11]  can  become  significant.  The  second  advantage  is  the  potential  to  measure 
the  nonlinear  behavior  of  the  tissue  itself.  This  behavior  is  itself  relevant  to  the 
physiological  functioning  of  several  tissues  (e.g.  arterial  walls,  cartilage),  and  may  be 
useful  in  distinguishing  potentially  cancerous  lesions.  It  is  the  latter  application  that 
directly  motivates  our  development. 

The  evidence  in  [12,13],  for  example,  indicates  that  the  degree  of  non-linearity  in  the 
stress-strain  relationship  of  soft  tissue  may  be  an  indicator  of  the  underlying  histology. 
In  [12],  the  mechanical  properties  of  soft  tissue  samples  of  the  breast  collected  from 
patients  during  surgery  were  studied  using  indentation  tests  at  various  strain  levels. 
The  authors  estimated  the  clastic  moduli  at  various  strain  levels  and  found  that  the 
cancerous  and  benign  breast  tissues  have  significant  differences  in  the  rate  of  increase 
of  stiffness  with  strain.  The  authors  of  [13]  tested  the  mechanical  properties  of  both 
breast  and  prostate  tissues  at  three  different  strain  rates  and  strain  levels.  They  also 
found  that  cancerous  tissues  are  much  stiffer  at  a  higher  strain  level  as  compared  to 
fat  or  normal  glandular  tissue.  Both  [12, 13]  make  their  measurements  at  times  scales 
associated  with  quasi-static  elasticity  imaging  (^  0.1  Hz  —  10 Hz),  and  report  negligible 
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viscoelastic  effects.  That  is,  the  observed  elastic  modulus  is  independent  of  loading 
frequency.  Of  these  various  applications,  our  primary  motivating  application  is  breast 
imaging,  with  the  aim  of  ultimately  improving  detection  and  differential  diagnosis  of 
breast  cancer. 

Relatively  few  attempts  have  been  made  to  account  for  non-linear  elastic  properties 
of  soft  tissue  in  material  property  reconstruction  [14-18].  Apart  from  the  exception 
of  [14]  only  constant  hyperelastic  parameters  are  recovered. 

In  [14]  the  authors  model  soft  tissue  as  a  linear  incompressible  elastic  material  accounting 
only  for  geometric  non-linearities.  The  relation  between  the  stress  and  strain  is  linear. 
The  authors  develop  an  inversion  equation  for  the  shear  modulus  in  an  integral  form 
and  solve  it  to  recover  shear  modulus  distributions.  The  displacement  fields  and  strain 
components  are  calculated  from  ultrasound  images  of  a  gel  based  phantom.  Accuarate 
results  were  obtained  for  recovered  shear  modulus  distributions  of  these  phantoms. 

In  [15]  the  authors  simulated  the  deformation  of  an  agar  gelatin  phantom  using  a 
commercial  finite  element  code  and  created  synthetic  RF  images  of  the  deformation 
of  the  phantom.  Using  displacement  data  calcualted  from  these  images,  the  relative 
strain  ratio  images  were  created.  Nonlinear  effects  were  observed.  That  is,  at  low  strain 
rates,  strain  contrast  was  observed  between  gelatin  (stiffer)  and  agar  (softer).  With 
increased  deformation,  the  agar  stiffened,  and  became  roughly  equal  to  the  stiffness  of 
gelatin.  This  caused  a  decrease  in  the  relative  strain  ratio.  With  further  increase  in  the 
load,  greater  differentiation  is  seen  between  the  agar  and  gelatin.  These  results  point  to 
the  fact  that  non-linear  properties  can  indeed  be  used  to  increase  differentiation  in  soft 
tissues.  Similar  studies  were  carried  out  by  the  authors  in  [17]. 

In  [16]  the  authors  develop  a  reconstruction  procedure  to  calculate  constant  parameters 
(or  their  combinations)  in  hyperelastic  material  modls  from  measured  force  displacement 
curves.  The  study  is  performed  both  for  phantoms  and  simulations.  A  similar  approach 
to  determine  the  hyperelastic  parameters  of  breast  tissue  samples  was  developed  in  [19] 
and  in  [18]  for  fitting  hyperelastic  models  with  constant  material  parameters  to  Treolar 
and  Treolar  and  Jones  experiments  on  natural  rubber. 

The  method  presented  below  was  developed  and  motivated  by  applications  in  breast 
imaging,  with  the  aim  of  ultimately  improving  detection  and  differential  diagnosis  of 
breast  cancer  [20].  In  these  applications  the  tissue  is  deformed  quasi-statically  while  it 
is  being  imaged  with  ultrasound.  Image  processing  yields  an  estimate  of  the  deformation 
pointwise  with  the  tissue,  at  several  “instants”  of  time.  The  time  scales  of  the  applied 
deformations  justify  the  incompessibility  assumption  and  neglect  of  inertial  and  viscous 
effects  [12,  13].  The  deformation  is  typically  measured  in  loading  only,  and  so  the 
quasi-elastic  [21]  model  is  justified.  We  therefore  assume  a  hyperelastic  model  of  tissue 
deformation.  The  aim  is  to  recover,  given  the  deformation  field  everywhere  in  the 
solid,  the  spatial  distribution  (in  the  material  or  reference  configuration)  of  the  tissue 
parameters  that  describe  the  assumed  hyperelastic  model. 
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1.1.  Organization  of  this  paper 

This  paper  is  organized  as  follows.  We  first  describe  the  strong  form  of  the  forward 
hyperelastic  problem  in  section  2.  This  serves  to  introduce  our  modeling  assumptions, 
definitions,  and  field  equations.  We  then  present  the  weak  form  which  leads  us  to  the 
discretization  and  computational  solution  of  the  forward  problem.  The  formulation  of 
the  inverse  problem  follows  in  section  3.  We  formulate  this  as  a  nonlinear  optimization 
problem.  Essential  to  its  practical  solution  is  the  efficient  evaluation  of  the  gradient. 
We  use  the  adjoint  method  and  a  smart  iteration-initialization  strategy  to  accomplish 
this,  as  described  in  section  4.  We  demonstrate  the  method  on  several  examples  in 
section  5.  Among  the  several  details  discussed  are  the  choice  of  hyperelastic  model, 
the  choice  of  the  regularization  functional,  the  selection  of  the  regularization  constant, 
the  effect  of  the  magnitude  of  the  applied  deformation  on  the  ability  to  recover  the 
nonlinear  parameters,  and  the  limits  and  appropriate  interpretation  of  linear  clastic 
reconstructions  in  hyperelastic  contexts. 

2.  The  Forward  problem 

2. 1 .  Strong  form 

In  this  section  we  describe  the  strong  form  of  the  equations  of  equilibrium  for  a  hy¬ 
perelastic  solid  undergoing  finite  deformations.  In  this  manuscript,  this  problem  is 
referred  to  as  the  forward  problem  of  hyperelasticity  or  simply  as  the  forward  problem. 
Given  appropriate  boundary  data  and  the  distribution  of  material  parameters  inside 
the  body,  these  equations  may  be  solved  to  determine  the  deformation  held  <f(X )  that 
maps  every  point  X  in  the  reference  or  material  configuration  to  a  point  x  in  the  cur¬ 
rent  or  spatial  configuration.  That  is  x  =  <p(X ).  The  displacement  held  is  given  by 
U(X)  =  <f>(X )  —  U(X).  This  is  illustrated  in  figure  (1)  For  a  complete  treatment  of 


Figure  1.  The  deformation  diagram 

this  problem  the  reader  is  referred  to  [22,23]. 

The  assumption  of  hyperelasticity  implies  that  there  exists  an  underlying  strain  energy 
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S 


^(t1,x2,x3;/3) 

dE 


(1) 


In  equation  (1),  E  is  the  Green-Lagrange  strain  tensor  and  S  is  the  second  Piola- 
Kirchhoff  stress  tensor.  The  Green-Lagrange  strain  tensor  may  be  expressed  as 
E  =  \{F7  F  —  1),  where  F  =  ' Vx(f)(X ),  is  the  gradient  of  the  deformation  <p  with 
respect  to  the  material  co-ordinates.  For  an  isotropic  solid,  the  strain  energy  density  0 
can  depend  on  the  deformation  gradient  F  only  through  the  list  of  invariants  (denoted 
by  X\ ,  Z2  and  Xf)  of  the  right  Cauchy-Green  strain  tensor  C  =  F1  F.  In  addition 
it  is  allowed  to  depend  on  a  vector  of  material  parameters  /3(X),  whose  values  may 
vary  spatially.  In  relation  to  soft  tissue  mechanics,  different  types  of  tissue  behavior  in 
response  to  applied  load  or  deformation  may  be  obtained  by  selecting  different  functional 
forms  of  the  strain  energy  function  0. 

The  equations  of  equilibrium  written  in  the  reference  configuration  Q0,  in  the 
absence  of  body  forces  are: 


Vx  •  (FS)  =  0  in  O0  (2) 

U  —  Q  on  T9  (3) 

FS  ■  N  =  XL  on  Tft  (4) 


In  equations  (2-4),  Ts  and  T/,  are  subsets  of  the  boundary  T  on  which  displacement  and 
traction  boundary  data,  Q  and  XL  respectively,  are  specified.  Note  that,  T  =  T/,  [jr9 
and  Thf^Tg  =  0,  that  is,  either  displacement  or  traction  boundary  conditions  must 
be  specified  on  the  entire  boundary.  In  equation  (4),  N  denotes  a  unit  normal  vector 
in  the  reference  configuration.  With  this  background  about  the  strong  form  of  the 
equilibrium  equations,  we  proceed  to  the  weak  form  of  the  hyperelasticity  equations 
and  their  computational  solution. 


2.2.  The  weak  form  and  its  numerical  solution 


We  begin  by  defining  functions  spaces  S  and  V,  for  admissible  solutions  and  weighting 
functions. 


S  =  {U  \Ui  e  H\n o);  Ui  =  Si  on  TJ  (5) 

V  =  {V\  Vi  EH1  (n0 );  Id:  =  o  on  TJ  (6) 


The  weak  formulation  of  the  equilibrium  equations  (2-4)  can  be  obtained  by  multiplying 
equation  (2)  by  a  weighting  function  W  G  V  and  integrating  by  parts.  Alternatively 
it  may  be  obtained  directly  by  minimizing  the  total  potential  energy  of  the  system.  In 
the  reference  configuration  it  is  given  by:  Find  a  displacement  U  G  S  such  that 


Wi  ■  Hi  dflo  =  o 


Th 


VW  G  V, 


A{W,Urf) 


(7) 
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A(W,U;/3)  =  [  WijFijSj 7 


dQn 


(8) 


The  weak  form  may  be  approximated  using  Galerkin’s  method  in  conjunction  with 
standard  finite  element  basis  functions  to  obtain  an  approximate  numerical  solution  of 
the  forward  problem.  This  solution  requires  the  linearization  of  the  semilinear  form  A, 
about  an  arbitrary  displacement  held  state  U  which  yields  a  bilinear  form  Ajj(-,  •;  f3,  U) 
given  by 


d 


Ajj(W ,  A U-  (3,U)  =  —  A(W,  U  +  eAC7 ;  (3) 


de^o 

/  Witj(5ikSjL  +  FuFkKCjjKL^U^L  dido- 
'do 


(9) 


In  equation  (9)  Cijkl  are  the  components  of  a  fourth  order  tensor  defined  as  follows: 


C-dS  ovC 

C~dE  C,JKL 


dS 


ij 


or  C 


dEKL  IJKL  dEudEKL 
Since  the  order  of  differentiation  does  not  matter,  it  is  clear  from  the  equation  above 
that  the  fourth  order  tensor  C  possesses  major  symmetry,  that  is  C-ukl  =  Cklij- 
Using  this  and  the  symmetry  of  the  stress  tensor  S  it  may  be  verified  from  (9)  that 
Ajj(-,  •;  /3,  U)  is  symmetric. 


d  20 


(10) 


2.3.  Numerical  solution 

The  solution  of  the  equations  of  equilibrium  for  a  solid  undergoing  finite  deformations 
is  significantly  harder  than  the  corresponding  linear,  infinitesimal  deformation  case. 
The  former  leads  to  a  set  of  nonlinear  algebraic  equations  which  may  not  converge  for 
large  applied  loads.  To  overcome  this  problem,  an  incremental  loading  strategy  is  often 
adopted  which  involves  applying  the  prescribed  boundary  conditions  in  several  (riioa(j) 
load  steps  and  using  nnewt(m  Newton  iterations  at  each  step.  This  leads  to  nioad  x  nnewton 
stiffness  matrix  assemblies  and  solves  as  compared  to  just  one  for  the  linear  case.  For 
our  applications,  we  have  found  that  nioad  ~  80  and  nnewton  ~  6,  which  makes  the 
solution  of  the  forward  problem  very  expensive.  Further,  since  the  formulation  of  the 
inverse  hyperelasticity  problem  we  propose  in  this  manuscript  requires  a  series  of  forward 
solves,  it  is  imperative  that  the  time  taken  to  solve  the  hyperelasticity  forward  problem 
be  kept  within  a  manageable  limit.  In  Section  (4.2),  we  describe  an  approach  which 
accomplishes  this. 


3.  The  formulation  of  the  inverse  problem 

We  formulate  the  inverse  problem  as  a  constrained  minimization  problem  and  solve  it 
using  a  quasi-Newton  optimization  algorithm.  Such  an  algorithm  requires  at  every 
iteration  the  gradient  of  the  objective  function  with  respect  to  the  optimization 
parameters.  The  gradient  of  the  objective  function  is  evaluated  efficiently  using  the 
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adjoint  of  the  linearized  hyperelasticity  equations.  Similar  approaches  for  solving  inverse 
problems  have  been  used  in  [4,24-26].  We  consider  a  generic  hyperelastic  material 
model,  which  depends  on  Np  material  parameters  (3  =  [/?i,  *  •'  • ,  /3np]t-  We  seek  the 
spatial  distribution  of  material  parameters  which  minimizes  the  difference  between  a  set 
of  n  measured  displacement  held  U™’1,  ■  ■  ■ ,  U"l,n  and  the  corresponding  set  of  predicted 
displacement  fields  U1,  The  predicted  displacement  helds  are  constrained  to 

satisfy  the  equations  of  equilibrium.  The  objective  function  is  the  square  of  the  L'2  norm 
of  the  difference  between  the  measured  and  predicted  displacement  helds.  We  have 
chosen  the  L 2  norm  as  it  avoids  the  need  to  differentiate  a  noisy  measured  displacement 
held  and  yields  a  maximum  likelihood  estimate  of  the  material  parameters  if  the  noise 
in  the  measured  displacement  helds  is  Gaussian  [27,  pp657-658]. 

The  statement  of  the  inverse  problem  is  :  Find  the  spatial  distribution  of  the 
material  parameters  / 3  such  that  the  objective  function  tt  (defined  in  equation  (11)) 
is  minimized  to  the  subject  to  the  constraint  that  each  deformation  held  satisfy  the 
equilibrium  equations,  that  is  (1-4)  or  (7)  and  (1). 

The  objective  function  n  is  given  by 

n  Np 

*  =  l'EwrmiF)-T(u^mi+1-Y,aM\\l  (u) 

p= i  j= i 

In  the  above  equation,  for  the  pth  measurement,  Um"p  is  the  measured  displacement 
held,  Up  is  the  corresponding  predicted  displacement  held,  T  is  a  tensor  which  selects 
the  appropriate  displacement  components  and  wp  represents  a  weighting  factor.  In 
addition,  for  each  material  property  held  /3j ,  a3  is  a  regularization  parameter  and  ||  •  H& 
is  the  regularization  norm.  The  significance  of  these  quantities  is  discussed  below. 

Remarks 

(i)  The  weighting  factors  wp  are  chosen  so  that  all  measurements  contribute  equally  to 
the  objective  function  tt.  This  ensures  that  there  is  no  bias  towards  measurements 
with  large  deformation.  For  example,  if  two  measurements  were  available  and  one 
of  these,  say  the  second  contained  much  larger  displacements.  Then  the  objective 
function  will  be  more  sensitive  to  the  second  held,  and  the  minimization  algorithm 
would  be  biased  in  favor  of  minimizing  this  held.  In  order  to  remove  this  bias,  it 
may  be  useful  to  select  the  weights  to  be  inversely  proportional  to  a  measure  of 
magnitude  of  the  displacement  helds.  That  is 

W\  =  (— )2  and  w2  =  1  (12) 

Qi 

where  q\ ,  q2  represent  the  magnitude  of  the  measured  displacements. 

(ii)  Often  only  one  component  of  any  measured  displacement  vector  Um,p  is  reliably 
known.  This  is  typically  the  case  when  the  imaging  system  has  better  resolution  in 
one  direction  than  the  other.  For  example,  when  ultrasound  is  used  to  measure 
the  deformation,  reliable  estimates  of  the  displacement  held  are  obtained  in  a 
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direction  parallel  to  the  axis  of  the  transducer.  In  such  cases  it  might  be  useful  to 
compare  only  the  axial  component  of  the  measured  displacement  with  the  predicted 
displacement.  In  our  formulation,  this  is  accomplished  by  the  tensor  T.  Assuming 
the  axial  component  is  along  the  unit  vector  e2,  an  appropriate  choice  for  T  is 

T  =  e2®  e2.  (13) 

(iii)  In  equation  (11)  oq  is  a  regularization  factor,  which  controls  the  trade  off  between 
matching  the  predicted  and  measured  displacement  fields  and  certain  desirable 
physical  characteristics  of  the  spatial  distributions  of  the  material  parameters 
(such  as  smoothness).  These  characteristics  are  selected  by  the  choice  of  the 
regularization  norm  ||  •  ||^.  The  choice  of  L 2  norm,  defined  by  (14)  ensures  that  the 
material  parameter  distribution  does  not  attain  unphysically  large  values.  The  H 1 
semi-norm  defined  in  equation  (15)  penalizes  oscillations  in  the  solution,  ensuring 
smoothness  in  the  recovered  material  parameter  distributions.  The  total  variation 
diminishing  (TVD)  regularization  norm,  defined  in  equation  (16),  does  just  that, 
but  does  not  penalize  jumps  in  the  solution.  In  this  norm,  c  is  a  small  parameter 
which  ensures  differentiability  of  the  norm  at  |  V/3®|  =  0.  In  this  manuscript  we 
have  chosen  c  =  0.1. 


mi 
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4.  Inverse  problem  solution  using  the  adjoint  method 


4.1.  Definition 


We  denote  the  variation  in  a  function  f(x)  in  the  direction  Sx  by  Sf  and  it  is  given  by: 


5f  =  DJ  ■6x=± 
de 


e^O 


f(x  +  eSx ) 


(17) 


Here  Dxf  is  the  directional  derivative  of  the  function  f(x). 


4-2.  Gradient  calculation  via  the  adjoint  method 

In  the  previous  section,  we  have  formulated  the  inverse  hyperelasticity  problem  as  a 
constrained  minimization  problem.  To  solve  this  problem  we  rely  on  a  quasi-Newton 
algorithm  which  requires  the  gradient  vector  at  each  iteration.  The  computation  of  the 
gradient  is  by  far  the  most  expensive  part  of  our  solution  algorithm.  In  this  section  we 
present  a  method  to  evaluate  the  gradient  efficiently. 

We  use  the  following  notation:  A  superscript  index  denotes  a  quantity  associated 
with  the  measurement  indicated  by  the  superscript.  For  example,  U'In'p  denotes  the  Ith 
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1:  Guess  a  homogeneous  distribution  for  each  of  the  material  parameters  /3i(X),  /32{X),  •  •  •  ,  div^lW) 

2:  repeat 

3:  Solve  forward  hyperelasticity  problems  corresponding  to  the  current  distribution  of  the  material 

parameters  using  the  previous  displacement  field  if  available 
4:  Solve  corresponding  linear  adjoint  problems 

5:  Evaluate  the  objective  function  and  its  gradient 

6:  Use  a  quasi-Newton  optimization  algorithm  to  update  the  guess  of  the  material  parameters 

7:  until  Either  the  objective  function  n  is  low  enough  or  the  maximum  number  of  iterations  have 
been  reached. 

Figure  2.  Solution  of  the  hyperelasticity  inverse  problem 


component  of  the  measured  displacement  field  corresponding  to  the  pth  measurement 
and  Up  denotes  the  Ith  component  of  the  predicted  displacement  held  corresponding  to 
the  same  measurement. 

We  note  that  in  the  objective  function  it,  (11),  each  U‘  depends  on  /3  through 
the  equilibrium  equations  (1)  and  (7).  We  account  for  this  constraint  using  Lagrange 
multipliers.  We  begin  by  writing  the  Lagrangian  C  corresponding  to  the  objective 
function  n  and  the  hyperelasticity  constraints. 

■  n  ..  Np 

CiW .  C/./3)  V  -£>, \\(T(Ut)  -  T(U™))\\1  +  - 

P=  1  j= 1 

n 

+  U‘\  0)  -  (W-,  W»)rj)  (18) 

i=p 

In  equation  (18),  Wp  is  interpreted  as  a  Lagrange  multiplier  enforcing  the  constraint 
of  hyperelasticity  for  the  pth  predicted  displacement  held.  In  the  following  development 
we  will  use  to  the  Lagrangian  to  derive  an  expression  for  the  derivative  of  the  objective 
function  with  respect  to  optimization  parameters. 

The  variations  in  the  Lagrangian  are  denoted  by  5C  and  are  given  by 

n  n 

=  DwpC  '  6WP  +  5Z  DUpC  '  6JJP  +  D(3C  '  ^  (19) 

p=  1  p=  1 

We  set  the  variations  in  C  with  respect  to  the  Lagrange  multipliers  Wp,  to  zero.  This 
yields  the  following  equation  for  each  Up  G  Sp, 

A(5WP,  Up ■  (3)  -  (5WP,  1-CP)  =  0  V5WP  G  Vp.  (20) 

Note  that  the  equation  above  is  identical  to  (7)  up  to  the  superscript  p.  That  is,  each 
Up  is  required  to  satisfy  the  forward  problem  for  that  particular  deformation  held. 
Equation  (20)  is  often  referred  to  as  the  primal  problem.  When  all  Up  are  selected  such 
that  equation  (20)  is  satished,  that  is  when  all  the  predicted  displacement  helds  Up 
satisfy  the  forward  elasticity  equations,  from  (1)  and  (18)  we  have, 


7T  =  C. 


(21) 
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Therefore  on  surfaces  satisfying  (20)  we  have 

5tt  =  SC.  (22) 

We  return  to  (19)  and  set  variations  in  C  with  respect  to  every  Up  equal  to  zero. 

That  is,  we  choose  values  for  the  Lagrange  multipliers  Wp  G  Vp  such  that 

DjjpC  ■  5UP  =  0  ydUp  G  Vp.  (23) 

Using  (18)  and  (8),  this  equation  yields 

Ajj(Wp,  8UP ;  A  Up)+  [  ( T(Up-Um’p))-T(5Up )  dU0  =  0  V 5UP  G  Vp(24) 

J 

where  the  bilinear  form  •;  A  Up)  is  given  by  (9).  Equation  (24)  which  is  used  to 

determine  Wp  is  the  adjoint  of  the  linearized  elasticity  equations.  Since  the  linearized 
operator  is  self-adjoint  (see  Section  2.2),  this  equation  may  be  written  as:  Find  Wp  G  Vp 
such  that 

Ajj(5Up,  Wp ;  /3,  Up)+  [  (T(Up-Um’p))-5Up  d<A0  =  0  V 5UP  G  Vp(25) 

J  Qo 

Note  that  when  Up  and  Wp  are  chosen  according  to  (20)  and  (25)  respectively,  the 
relation  (22)  holds  and  the  first  two  terms  in  (19)  vanish.  As  a  result  we  have 

8n  =  Dp£  •  5(3 


Equation  (26)  represents  the  change  in  n  due  to  an  infinitesimal  change  in  the  material 
parameters  in  the  direction  of  5/3. 

To  derive  an  expression  for  the  gradient  for  the  finite  dimensional  problem,  we 
expand  the  material  parameters  in  terms  of  a  basis, 

N 

ft  =  £  MA(X)bjA  (27) 

A=  1 

In  the  above  equation,  N  represents  the  number  of  basis  functions,  Ma(X),  used 
to  represent  each  material  parameter  and  bjA  is  the  coefficient  of  the  j'-tli  material 
parameter  corresponding  to  the  Atli  basis  function.  As  a  result  variations  in  the  material 
parameters  are  given  by 

N 

5(5j  =  YJMA(X)5b3A  (28) 

A= 1 

where  5b3A  represents  variations  in  bjA ■  Using  (28)  in  (26)  we  have 

Np  N 

5n  =  EE  fij  a  5b  jA 

j= 1  A=  1 


(29) 
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where  gjA  is  identified  as  the  gradient  of  n  with  respect  to  the  material  parameters  with 
respect  to  the  Ath  coefficient  for  the  jth  material  parameter.  From  (26)- (29)  we  conclude 

n  „ 

9,  a  =  £  /  WfrFtJ 

P=1  v 

Note  that  in  equation  (30)  there  is  no  sum  on  the  index  j .  Thus  in  order  to 
evaluate  gjA  at  every  iteration,  we  first  solve  n  non-linear  forward  problems  (20)  to 
determine  Up  p  =  1,2 ,■■■  ,n,  then  solve  n  linear  adjoint  problems  (25)  to  determine 
Wp  p  =  1, 2,  •  •  • ,  n  and  use  this  in  (30).  Our  algorithm  to  solve  the  inverse  problem  is 
described  in  Table  (2). 

Remarks 

(i)  Increasing  computational  speed  of  the  inverse  problem:  The  solution  time 
of  the  inverse  problem  with  the  adjoint  method  is  dominated  by  the  computational 
cost  of  solving  the  nonlinear  forward  problem  (20)  at  each  iteration.  The  high 
computational  cost  of  the  forward  problem  arises  because  it  needs  to  be  solved 
using  small  load  increments  and  performing  several  Newton  iterations  for  each 
increment.  To  speed  up  the  solution  of  this  problem,  and  hence  the  computation 
of  the  gradient,  we  perform  the  load  increments  only  the  first  time  the  forward 
problem  is  solved  (the  first  iteration  of  the  optimization  algorithm).  For  subsequent 
iterations,  we  apply  the  entire  load  in  one  iteration.  This  requires  a  very  good  guess 
of  the  displacement  field  everywhere  in  the  domain,  which  is  the  displacement  field 
corresponding  to  the  previous  iteration  of  the  optimization  algorithm.  We  have 
found  that  this  field  is  close  enough  to  ensure  quick  convergence.  Thus,  after  the 
first  iteration,  the  forward  elasticity  problem  is  solved  very  efficiently  (typically  in 
about  6  Newton  iterations). 

(ii)  We  have  not  addressed  the  uniqueness  of  the  recovered  spatial  distribution  of 
material  parameters.  Indeed,  there  could  be  more  than  one  set  of  material  property 
distributions,  perhaps  even  entire  families  of  distributions,  that  could  minimize 
the  objective  function  (11).  For  the  inverse  problem  for  the  linear,  infinitesimal 
deformation  elasticity  case,  it  has  been  demonstrated  that  the  uniqueness  of 
the  solution  depends  on  the  boundary  conditions,  and  that  for  problems  with 
prescribed  displacement  boundary  condition  a  single  measured  displacement  field 
is  not  sufficient  to  ensure  uniqueness  [28,29].  It  has  also  been  shown  that  multiple 
measured  deformation  fields  (for  example  deforming  the  same  elastic  specimen  in 
compression  and  in  shear)  can  be  used  to  restrict  the  set  of  solutions  for  the  spatial 
distributions  of  material  parameters  and  render  the  solution  unique.  We  believe 
that  a  similar  analysis  can  be  performed  in  the  finite  deformation,  hyperelastic  case 
for  a  specific  form  of  the  strain  energy  density  function  or  to  identify  the  coefficients 
oq(cc),  a2(x)  in  the  general  form  of  the  Cauchy  stress  tensor  and  the  equilibrium 
equation  in  the  current  configuration  as  below: 

T  =  -pi  +  a1{x)B  +  a2{x)B^1  (31) 


D 


Ma{X)  dn o  +  a^,  MA{X))b.  (30) 
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Wx  T  =  0  (32) 

It  is  also  clear  that  for  non-linear  behavior  of  tissue  distinct  displacement  fields 
can  also  be  obtained  by  using  different  levels  of  compression  of  the  elastic  medium. 
We  are  in  the  process  of  extending  the  work  in  [28,  29]  to  the  finite  deformation, 
hyperelastic  case. 

(iii)  The  strong  form  of  the  equations  of  equilibrium  (2)  can  be  treated  as  a  set  of 
equations  for  the  unknown  material  parameters  (3  with  the  measured  displacement 
held  U  as  a  known  quantity.  These  equations  may  then  be  solved  for  the  material 
parameters  /3.  This  would  be  the  so-called  “direct  approach”  to  solving  the  inverse 
elasticity  problem.  The  advantage  of  this  approach  is  its  computational  ease 
(solving  a  single  nonlinear  partial  differential  equation).  However,  this  approach 
also  involves  computing  the  second  derivative  of  the  displacement  held,  which 
may  be  inaccurate  in  the  presence  of  noise.  In  addition,  it  relies  on  knowing  all 
components  of  the  measured  displacement  held,  which  may  not  be  available  in 
practice. 

(iv)  The  adjoint  method  of  gradient  calculation  requires  only  one  solve  of  the  forward 
hyperelasticity  problem  and  one  solve  of  the  adjoint  problem  to  compute  the 
gradient  of  the  objective  function.  This  means  that  the  cost  to  compute  the  gradient 
of  the  objective  function  is  essentially  independent  of  both  the  number  of  material 
parameters  in  the  hyperelastic  constitutive  relation  and  the  number  of  parameters 
used  to  discretely  represent  each  parameter. 

(v)  In  this  work,  we  have  not  considered  the  question  of  possible  loss  of  material  stability 
and  its  effect  on  the  inversion  algorithm.  The  inequalities  of  Baker-Ericksen, 
Coleman-Noll,  and  strong  cllipticity  need  to  be  satisfied  for  each  iteration  in  the 
inversion  algorithm.  These  inequalities  are  motivated  by  physical  considerations 
such  as  “stress-increases  with  strain”.  These  inequalities  may  imply  certain 
additional  inequalities  among  the  material  parameters  appearing  in  the  hyperelastic 
relation  which  may  change  with  the  deformation.  Not  satisfying  these  inequalities 
may  lead  to  problems  of  convergence  for  the  finite  element  method  and  may  lead  to 
unphysical  parameters  in  the  constitutive  relation  .  The  reader  is  referred  to  [23] 
for  a  discussion  of  these  inequalities. 

5.  Numerical  examples 

In  this  section  we  present  two  dimensional  examples  in  which  we  reconstruct  the  of  dis¬ 
tribution  non-linear  parameters  from  synthetic  displacement  fields  with  noise.  We  first 
assume  a  distribution  of  material  properties  corresponding  to  a  stiff  circular  inclusion 
in  a  homogeneous  background.  We  then  solve  two  forward  hyperelasticity  problems 
corresponding  to  this  distribution  of  material  properties  and  obtain  the  corresponding 
displacement  fields.  In  the  examples  considered,  one  of  the  displacement  fields  is  with 
large  applied  strain  (20%)  and  the  other  displacement  held  is  with  small  applied  strain 
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(1%).  The  Veronda-Westman  constitutive  relation  is  used  to  describe  the  constitutive 
behavior  of  soft  tissue.  This  strain  energy  function  and  associated  quantities  for  this 
relation  are  presented  in  equations  (33,34,35).  Additive  white  Gaussian  noise  is  then 
added  to  the  synthetic  displacement  fields  thus  generated  using  the  AWGN  function 
in  Matlab,  to  simulate  experimentally  acquired  data.  Such  a  noise  level  corresponds 
to  40dB  level  noise  in  the  measured  displacement  field  and  is  believed  to  be  realistic 
representation  of  the  level  of  noise  in  an  ultrasound  system  [6]. 

We  then  recover  the  material  properties  using  only  one  component  (the  axial  compo¬ 
nent)  of  the  displacement  fields.  We  start  with  a  homogeneous  guess  of  the  material 
properties  and  then  improve  on  this  homogeneous  guess  by  using  a  quasi  Newton  opti¬ 
mization  approach.  The  gradient  of  the  objective  function  is  calculated  efficiently  using 
the  adjoint  method  as  described  in  section  (4.2).  For  this  purpose,  we  have  developed  a 
hyperelastic  finite  element  based  code  for  the  calculation  of  the  hyperelastic  parameters. 
The  results  presented  were  obtained  on  a  1.7Hz  Pentium  Xeon  machine. 

We  begin  by  describing  the  constitutive  relation  used  to  describe  the  mechanical  behav¬ 
ior  of  soft  tissue  in  section  (5.1)  and  discuss  the  significance  of  the  material  parameters 
appearing  in  it.  We  then  linearize  the  material  tangent  modulus  and  show  that  in  the 
limit  of  small  deformations  the  deformation  is  affected  only  by  the  parameter  p,  and  this 
then  motivates  the  use  of  a  small  strain  deformation  field  (1%)  in  the  reconstruction. 
We  then  specialize  the  three  dimensional  constitutive  law  to  two  dimensions  assuming 
plane  strain  deformations  in  section  (5.1.2).  We  then  describe  the  boundary  conditions 
and  the  geometry  of  the  specimen  which  is  deformed  in  section  (5.2).  We  close  with 
figures  showing  ideal  and  recovered  distributions  of  the  material  parameters  in  section 
(5.3). 

5.1.  Constitutive  relation  and  significance  of  the  parameters 

In  the  non-linear  elasticity  regime,  the  constitutive  relation  is  of  paramount  importance 
because  this  is  the  relation  that  governs  the  behavior  of  tissue  in  response  to  applied 
deformation  or  load.  When  the  tissue  behavior  is  assumed  to  be  hyperelastic,  then 
different  forms  of  tissue  behavior  and  thus  different  tissue  models  can  be  chosen  by 
selecting  different  forms  of  the  stored  energy  function  0(X1,X2,X3,  (3).  Different  strain 
energy  functions  have  been  used  in  the  literature  relating  to  elasticity  imaging  to  describe 
tissue  deformation.  For  example,  to  characterize  the  behavior  of  soft  tissue,  the  Arruda- 
Boyce  relation  [30]  has  been  used  by  Li  et  al  in  [31]  to  model  the  behavior  of  soft  tissue 
in  an  indentation  experiment.  The  Mooney- Rivlin  strain  energy  function  has  been  used 
in  [32]  to  calculate  hyperelastic  properties  of  soft  tissues  using  a  truncated  expansion 
in  the  strain  invariants.  The  same  form  has  been  used  in  [16].  In  this  paper,  we  have 
selected  the  Veronda-Westmann  relation  [33]  characterized  by  equations  (33,34,35),  since 
soft  tissues  exhibit  approximately  exponential  stress  strain  response.  This  form  of  the 
strain  energy  function  was  used  by  Nienhuys  in  [34]  to  model  the  cutting  of  soft  tissue. 
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In  this  work  we  recover  the  parameters  7  and  p.  The  parameter  7  is  the  non-linearity 
in  the  relation  and  the  parameter  p  corresponds  to  the  shear  modulus  at  zero  strain.  (3 
is  a  large  parameter  that  penalizes  compressible  deformations. 


=  —  (e7^1-3)  -  1)  -  ^(X2  -  3)  +  /3(X3  -  1  -  ln(l3)) 

7  2 

(33) 

Su 

=  2  +  £(C„ 

+  C/;)  +  ^(^3  — +  0//) 

(34) 

Cijkl 

=  4  +  ^(SIKSJL  +  <W„) 

(35) 

+  f  feeder?  -  (i3  -  i)(c„;c;f + cj^cy)' 

5.1.1.  Linearization  of  the  material  tangent  modulus  For  small  deformations,  the  ma¬ 
terial  tangent  modulus  can  be  linearized  about  the  zero  strain  state  as  follows. 

CljKL  =  4(7^t  —  +  (3)5kl8iJ  +  p(8ik8jl  +  SjrSli )  (36) 

In  equation  (36),  /3  is  a  large  parameter  f3  7/x  —  Thus,  the  deformation  at  small 
strains  is  effectively  controlled  by  variations  in  coefficient  (p).  in  front  of  the  second  term 
(5ik5jl  +  bjicbu)-  In  solving  the  inverse  problem,  we  can  take  advantage  of  this  fact 
by  recovering  the  parameter  p  from  a  small  strain  deformation  field.  The  distribution 
of  the  parameter  7  is  recovered  by  incorporating  a  large  strain  deformation  field  in  the 
objective  function. 

5.1.2.  Specialization  to  two  dimensions  For  computational  purposes,  we  specialize  the 
three  dimensional  Veronda-Westman  constitutive  relation  (33,34,35)  to  two  dimensions 
by  assuming  plane  strain  deformations.  Then  form  of  the  Veronda-Westmann  relation 
changes  to 

<t>  _  —  (eT(Xl~2)  -  1)  -  ^(X2  +  -  3)  +  /3(J2  -  1  -  ln(X2))  (37) 

7  2 

Here  X!,X2  are  the  invariants  of  the  two  dimensional  right  Cauchy-Green  strain  tensor 

c  =  ftf. 


S/j 

=  2  [(^Ii-2)-|x1-|)<y„  +  ^(c/J 

+  Cj/)  +  ^(X2  -  1)(G/JT  +  Cj/) 

(38) 

Cijkl 

=  4[7/xe7(Xl-2)^L(57J  -  |W/j  +  ^(Wjl  +  <Wl/) 

(39) 

+  2  2 X2Ca-^C,jjt  —  (X2  —  1)(C'/^C'JJ  +  CjkCjT) 
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5.2.  The  problem  configuration,  the  deformation  fields  used:  1%  and  20%  strain 

We  have  used  two  deformation  fields  to  recover  the  distributions  of  parameters  7,  p  in 
the  Veronda-Westman  constitutive  relation.  The  deformation  field  chosen  are  1)  at  1% 
compression  2)  20%  compression.  These  deformation  fields  are  chosen  because  at  small 
deformations  only  the  parameter  p  which  is  the  shear  modulus  at  zero  strain  state  affects 
the  deformation  of  the  material.  At  larger  strains  the  parameters  7  and  p  both  affect  the 
deformation.  Thus  both  the  parameters  can  be  recovered  from  the  deformation  fields. 
In  both  deformations,  the  specimens  are  deformed  according  to  the  schematic  diagram  3. 


5.3.  Two  dimensional  examples:  figures  of  parameter  recovery 

We  present  the  following  results  for  the  example  problem  show  in  figure  (3). 

(i)  Ideal  distributions  of  p  and  7. 

(ii)  Recovered  distributions  of  p  and  7. 

(iii)  Line  plots  of  the  material  properties  through  sections  AA  and  BB  shown  in  figure 

(4). 


Problem  Size  (elements) 

Number  of  iterations 

Solution  time 

30X30 

100 

~  2700  seconds 

5-4-  Problem  description 

Figure  (3)  shows  a  schematic  of  the  problem  solved  in  the  reference  configuration. 
The  nonlinearly  elastic  medium  obeying  the  Veronda  Westman  relation  is  described 
in  equations  (33,34,35),  is  a  rectangular  domain  with  its  top  edge  held  fixed  in  the 
y-direction  (free  to  move  in  the  x  dirction)  and  displacement  boundary  conditions 
(compression)  applied  to  the  opposite  edge  in  the  y-dirction  with  the  nodes  being  free 
to  move  in  the  x-direction.  The  sides  are  traction  free.  Two  levels  of  compression  are 
applied  to  the  boundary  (1%  and  20%),  to  yield  two  different  displacement  fields.  The 
time  taken  to  generate  these  two  deformation  fields  is  30  minutes  or  15  minutes  each. 
That  is,  the  time  taken  for  one  complete  forward  solve  is  15  minutes. 

5.Ihl.  Ideal  ( Target )  distributions  of  p  and  7  Figure  (4)  shows  the  distributions  of  the 
material  parameters  p  and  7  that  where  used  to  generate  synthetic  displacement  data. 
Thus,  they  represent  the  ideal  spatial  distributions  of  material  parameters  that  can  be 
recovered  from  the  displacement  data. 

5.4-2.  Parameter  recovery  at  1%  noise  In  figures  (6,7,8)  we  show  the  recovered 
distributions  of  p  and  7  with  various  regularization  norms. 
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(a)  Target  distribution  of  7 


Figure  3.  Example  problem  used  for  numerical  simulations 


(a)  Target  distribution  of  7  (b)  Target  distribution  of  y 


Figure  4.  Target  distribution  of  7  and  y 


(a)  Stress  strain  curve  in  compres-  (b)  Incremental  elastic  modulus  in 
sion  compression 

Figure  5.  Response  Veronda  Westman  material  in  compression 


5. 5.  Lineplots  of  the  material  properties  through  a  section  in  their  center 

Figure  (9)  shows  variations  of  material  properties  along  the  vertical  sections  AA  and  BB 
which  are  shown  in  (4).  The  plots  how  that  in  this  case,  to  recover  /i,  the  L 2  norm  and 
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(a)  Recovered  distribution  of  7  (b)  Recovered  distribution  of  y 

Figure  6.  Recovered  distributions  of  7  and  y  with  the  TVD  norm 


5  10  15  20  25  30  5  10  15  20  25  30 


(a)  Recovered  distribution  of  7  (b)  Recovered  distribution  of  y 

Figure  7.  Recovered  distributions  of  7  and  y  with  the  HI  semi  norm 


(a)  Recovered  distribution  of  7  (b)  Recovered  distribution  of  y 

Figure  8.  Recovered  distributions  of  7  and  y  with  the  L2  norm 


the  H1  semi  norm  are  the  most  effective  norms  in  recovering  the  parameter  contrast. 
The  recovery  with  the  L 2  norm  is  non-smooth  and  the  inclusion  is  of  the  wrong  size.  The 
TVD  norm  recovers  roughly  the  same  contrast  ratio  albeit  with  different  average  and 
peak  values.  This  can  be  seen  in  Figure  (6)  in  which  the  entire  background  is  elevated 
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above  the  starting  guess.  Figure  (9)  also  shows  similar  plots  for  7.  In  this  case  the 
performance  of  the  TVD  norm  is  seen  to  be  good  in  recovering  the  sharp  edges.  This  is 
not  surprising  since  the  underlying  parameter  distributions  have  sharp  edges,  and  the 
TV D  norm  is  expected  to  yield  sharp  edges  in  reconstructions.  The  parameter  recovery 
for  7  with  the  L 2  norm  is  comparable  in  contrast  to  the  TVD  norm  but  the  parameter 
recovery  is  not  smooth.  The  H 1  —  semi  norm  yields  a  smooth  result  for  the  parameters. 

6.  Conclusions  and  future  work 

In  this  work  we  have  developed  an  efficient  iterative  formulation  for  the  non-linear 
elasticity  imaging  inverse  problem.  The  method  relies  on  the  adjoint  method  of 
optimization  to  calculate  the  gradient  of  a  cost-function  efficiently.  Representative 
numerical  examples  were  presented.  These  examples  show  that  it  is  feasible  to 
reconstruct  material  parameter  distributions  for  hyperelastic  materials  in  the  presence 
of  noise.  With  a  good  starting  displacement  guess  for  the  inverse  problem,  the  solution 
time  of  the  non-linear  elasticity  imaging  inverse  problem  using  the  adjoint  method  is 
comparable  to  the  solution  time  of  the  linear  elasticity  imaging  inverse  problem  using 
the  adjoint  method. 

It  must  be  emphasized  that  there  are  numerous  other  issues  involved  in  non-linear 
elasticity,  which  maybe  addressed  in  future  work.  These  include,  tests  on  experimental 
and  clinical  data,  the  importance  of  using  the  correct  constitutive  law  (or  the  effect  of 
using  an  incorrect  constitutive  law)  as  the  model  for  the  tissue  behavior,  feasibility  of 
recovering  non-linear  elastic  properties  when  small  strains  are  used  in  the  deformation 
process,  the  feasibility  of  using  different  tissue  constitutive  relations  in  different  parts 
of  the  domain,  possibly  adaptively,  or  the  feasibility  of  performing  reconstructions  for 
anisotropic  non-linear  elasticity. 
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(a)  Line  plots  representing  variation  of  7  along  a  vectical  section 
through  the  center 


(b)  Line  plots  representing  variation  of  y  along  a  vectical  section 
through  the  center 

Figure  9.  Line  plots  representing  variation  of  7  and  y  along  a  vertical  section  through 
the  center 
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Appendix  A.  Selection  of  the  regularization  parameters  an  and  a?2 

In  this  section,  we  describe  in  detail  the  selection  of  the  regularization  parameters  cq 
and  0:2  appearing  in  equation  ??.  Some  regularization  is  necessary  to  be  added  to  the 
problem  because  of  the  presence  of  noise  in  the  experimentally  acquired  data  or  to  ensure 
uniqueness  of  solution.  The  regularization  parameters  oq  and  02  are  calculated  for  the 
two  deformation  fields  at  1%  and  at  20%  according  to  Morozov’s  discrepancy  principle. 
The  reader  is  referred  to  [37]  for  a  review  of  this  and  other  regularization  techniques 
such  as  the  L-curve  method.  The  treatment  of  regularization  parameter  selection  closely 
follows  our  treatment  of  the  problem  in  [4],  We  begin  by  defining  the  noise  A  in  the 
measurement: 

A  =  || T(Umi)  -  T(£/I)||/||T(t/mi)||  (A.l) 

Here,  T(Umi )  represents  the  measured  component  of  the  ith  displacement  field.  U 
represents  the  noiseless  measurement.  Thus,  || T(Umi)  —  T(Ul)\\  is  the  deviation  from 
the  perfect  measurement.  Then,  cq  should  be  chosen  to  be  the  largest  real  number  that 
allows: 

|| T(Umi  -  IP)  ||  =  C\\T(Umi  -  U"d)  ||  (A. 2) 

In  the  above  equations,  C  is  a  real  number  which  is  chosen  such  that  C  ~  1.  This  means, 
that  the  difference  between  the  measured  and  predicted  displacement  field  should  be 
roughly  the  same  as  the  difference  between  the  measured  displacement  held  and  the 
noiseless  displacement  held.  In  other  words,  one  should  not  match  the  measured  and 
displacement  helds  to  a  greater  extent  than  the  level  of  noise  present  in  the  problem. 
Choosing  a  large  value  of  cq  increases  the  importance  of  the  regularization  term  in 
relation  to  \\T{Uim  —  Ul)  ||  and  gives  a  solution  for  the  material  parameter  distribution 
that  is  dominated  by  the  characteristics  chosen  by  the  choice  of  the  regularization  norm. 
However,  this  may  mean  that  the  measured  and  predicted  displacement  helds  are  not 
well  matched.  Lowering  the  regularization  parameter  a  reduces  \\T(Umi  —  Ul)  ||. 
Another  way  of  looking  at  Morozov’s  discrepancy  principle  is  the  following.  We  use  the 
triangle  inequality  to  write: 

||T(LT)  -  T{IT) ||  <  ||T(LT)  -  T(Umi)\\  +  \\T{Umi)  -  T(IT)\\  (A.3) 

In  the  above  equation,  the  left  hand  side  represents  the  total  error  in  the  reconstruction 
if  perfect  data  is  available.  The  right  hand  side  is  the  sum  of  the  error  in  reconstruction 
||T(LT)  —  T{Umi)  ||  and  error  in  measurement  \\T(Umi)  —  T(Ul)  ||  for  the  ith  deformation 
held.  Morozov’s  discrepancy  principle  states  that  the  regularization  parameter  should 
be  chosen  so  that  a  balance  between  the  reconstruction  error  and  the  measurement  error 
is  achieved. 

We  now  describe  how  cq  is  calculated  for  a  known  noise  level  A.  Substituting  equation 
(A.l)  in  (A. 2)  we  obtain  the  following  equation 

|| T{Umi  -  Ul)\\  =  CA\\T(Umi)\\ 


(A.4) 
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Quantity 

Field  1 

Field  2 

\\T(Um-U)  || 

3.2411E-3 

7.061E-3 

\\T(Um)\\ 

0.5752 

11.38 

A 

0.01 

0.01 

C 

0.5634 

0.6205 

Table  Bl.  Calculation  of  the  constant  C  (see  equation  A. 4)  for  the  TYD  norm 


Quantity 

Field  1 

Field  2 

\\T(Um  —  U)\\ 

3.36E-3 

0.10655 

\\T(Um)\\ 

0.5752 

11.38 

A 

0.01 

0.01 

C 

0.6382 

0.9363 

Table  B2.  Calculation  of  the  constant  C  (see  equation  A. 4)  for  the  HI  semi  norm 


Quantity 

Field  1 

Field  2 

\\T(Um  —  U)\\ 

4.42E-3 

0.0999 

\\T(Um)\\ 

0.575 

0.13211 

A 

0.01 

0.01 

C 

0.7729 

1.16 

Table  B3.  Calculation  of  the  constant  C  (see  equation  A. 4)  for  the  L2  norm 


The  term  A||T(C/m*)||  can  be  easily  calculated  once  A  is  known.  The  term  1 1 T ( Urni  — 
U1  )||  can  be  calculated  assuming  a  value  for  a  and  solving  the  hyperelasticity  inverse 
problem.  With  this  knowledge,  C  can  be  calculated.  We  then  repeat  this  process, 
adjusting  the  values  of  a  until  C  ~  1. 

Appendix  B.  Selection  of  the  regularization  parameters  for  this  example 

The  regularization  parameters  in  this  example  were  adjusted  so  that  C  ~  1.  This 
calculation  is  shown  in  tables  (B1,B3) 
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1.  Introduction 

Elastography  refers  to  a  collection  of  imaging  techniques  that  allow  mechanical  strain 
distributions  to  be  imaged  and  noninvasively  quantified  in  vivo. 

The  time  scales  over  which  the  tissue  response  is  typically  measured  ranges  from 
about  a  millisecond  (the  typical  duration  for  a  radiation  force  “push  pulse”  to  about 
one  second,  (the  typical  time  scale  of  freehand  quasistatic  compression  used  in  strain 
imaging  [1,  2]).  Magnetic  resonance  elastography  and  sonoelasticity  imaging  typically 
use  time-harmonic  excitations  with  frequencies  in  the  range  of  102Hz. 

The  mechanical  responses  of  soft  tissues  and  tissue  mimicking  gels  observed  under 
transient  excitations,  be  they  radiation  force  or  time-harmonic  excitations,  show 
a  predominantly  elastic  component  as  well  as  a  small  viscoelastic  component.  In 
quasistatic  deformations,  on  the  other  hand,  the  strain  fields  are  typically  observed 
for  about  a  second  and  are  interpreted  within  the  context  of  linear  (or  rarely  nonlinear) 
elasticity.  That  is,  the  tissue  response  is  assumed  and  observed  to  be  approximately 
purely  elastic. 

Soft  tissue  is  widely  recognized  as  having  both  fluid  and  solid  phases  which 
can  move  independently  of  each  other.  Furthermore,  the  fluid  exists  within  several 
“compartments”  of  the  soft  tissue,  notably,  the  vasculature  (including  both  the  hemal 
and  lymphatic  vessels)  and  the  extravascular  space.  Of  course,  due  to  permeability 
of  the  microvessels  in  both  vascular  networks,  fluid  is  often  exchanged  between  these 
compartments.  It  is  recognized  that  fluid  flow  leads  to  a  stress  relaxation  at  fixed 
strain  (or  conversely,  a  strain  relaxation  at  fixed  stress).  It  is  reasonable  to  conjecture 
then,  that  by  measuring  the  spatio-temporal  patterns  of  strain  in  a  strain-relaxation 
type  of  experiment,  the  effects  of  fluid  flow  can  be  visualized  and  measured.  Indeed, 
recent  experiments  on  a  poroelastic  tissue  mimicking  phantom  have  demonstrated  the 
ability  to  image  the  effects  of  fluid  flow  on  spatio-temporal  strain  patterns,  and  to 
interpret  those  effects  within  the  biphasic  [3]  or  Biot  poroelasticity  theory. 

The  linear  “biphasic  theory”  [3,  4]  can  be  regarded  as  a  special  case  of  Biot 
poroelasticity;  the  special  case  being  that  of  having  two  incompressible  phases.  It 
has  been  very  successful  at  modeling  the  fluid-elastic  coupling  in  cartilage  [4,  3]. 
Cartilage  tends  to  be  avascular,  however,  and  so  fluid  resides  only  in  the  “extravascular 
compartment” . 

A  different  model  for  the  mechanics  of  vascularized  soft  tissue,  which  includes  the 
effects  of  fluid  flow  and  the  possibility  of  exchange  between  fluid  compartments  was 
proposed  by  R.  Skalak,  RK  Jain,  and  coworkers.  The  model  was  originally  developed 
in  a  rather  general  context  to  capture  effects  of  fluid-elastic  coupling  in  soft  tissues, 
but  was  then  applied  to  describe  perfusion  and  drug  delivery  in  solid  tumors.  It  has 
since  been  applied  in  [5],  and  validated  in  an  experimental  model  in  [6]. 

Our  motivation  for  this  work  stems  from  the  question  Can  techniques  from 
elastography  be  used  to  image  and  quantify  interstitial  fluid  flow  in  soft  tissues  from 
spatio-temporal  patterns  of  elastic  strain?  To  answer  this,  we  use  the  mathematical 
model  of  [5]  in  conjunction  with  finite  element  modeling  to  predict  the  effects  of  fluid 
flow  on  the  spatio-temporal  patterns  of  soft-tissue  elastic  strain  under  a  variety  of 
physiological  conditions.  The  magnitude  of  the  strain  effects  and  their  time  scales 
dictate  the  measurability  of  the  effects  of  fluid  flow.  Simulations  relevant  to  a 
quasistatic  elasticity  imaging  for  the  characterization  of  fluid  flow  in  solid  tumors  are 
emphasized  here.  In  this  context,  the  following  questions  are  specifically  addressed: 

(i)  How  do  characteristics  of  tumor  microvasculature  effect  strain  relaxation  in  solid 
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tumors? 

(ii)  How  does  the  spatio-temporal  strain  pattern  depend  on  the  relative  importance 
of  fluid  flow  within  the  extravascular  compartment  versus  fluid  exchange  between 
the  vascular  and  extravascular  compartments? 

(iii)  How  does  the  choice  of  boundary  conditions  effect  the  spatio-temporal  patterns 
of  strain? 

In  the  following  “Methods”  section,  we  describe  a  mathematical  model  that  is 
used  to  address  the  questions  enumerated  above.  We  present  an  exact  analytical 
solution  of  this  model  to  be  used  to  develop  intuition  and  to  serve  as  a  check  on 
a  finite  element  implementation.  We  then  describe  four  computational  experiments 
designed  to  answer  the  questions  raised  above.  This  is  followed  by  Results,  Discussion 
and  Conclusions.  In  the  Appendix,  we  include  a  derivation  of  the  field  equations  and 
the  exact  analytical  solution. 

2.  Methods 

2.1.  Mathematical  Model 

We  use  the  mathematical  model  described  in  [5],  and  derived  in  Appendix  A.  The 
model  treats  the  interstitial  space  as  a  biphasic  material,  and  incorporates  fluid 
exchange  between  the  interstitial  compartment  and  the  microvasculature.  It  is  this 
fluid  exchange  that  distinguishes  this  model  from  a  biphasic  model  for  nonvascular 
soft  tissues  such  as  cartilage. 

The  assumptions  that  go  into  the  model  are  small  strains,  small  vascular  space, 
Starling’s  law  for  (transient)  fluid  transport  across  the  vessel  wall,  Darcy’s  law  for  fluid 
flow  through  the  interstitial  compartment,  and  Hooke’s  law  for  the  elastic  response. 
We  further  assume  the  deformation  takes  place  slowly  enough  that  inertia  can  be 
neglected.  Under  these  conditions  (see  Appendix  A  for  details),  the  solid  displacement 
vector  u  and  interstitial  fluid  pressure  p  are  related  by: 

V  •  u  —  V  •  [«Vp]  +  XP  =0  (1) 

V  •  [-pi  +  AV  ul  +  2 pVSymu\  =  0  (2) 

Equation  (1)  represents  a  combination  of  the  conservation  of  fluid  mass  in  the 
interstitium,  with  the  momentum  equation  for  the  fluid  phase.  Equation  (2)  represents 
the  balance  of  total  linear  momentum  in  the  tissue.  The  symbols  that  appear  in 
equations  (1)  and  (2)  are  defined  as  follows:  V  is  the  gradient  operator;  I  is  the 
identity  tensor;  u  =  du/dt  is  the  solid  phase  velocity;  k  is  the  interstitial  permeability 
that  governs  the  ease  by  which  fluid  percolates  through  the  interstitium;  A  and  p 
are  the  drained  elastic  Lame  parameters  of  the  interstitium;  and  \  is  the  average 
microfiltration  coefficient,  given  by  %  =  xv  +  Xl >  with  xv  =  LpyV  and  Xl  =  LpySh  ; 
Lp  (resp.  Tpi)  is  the  hydraulic  conductivity  of  the  hemal  (resp.  lymphatic)  capillary 
wall,  Sy/V  (resp.  Sl/V )  is  the  surface  area  of  the  hemal  (resp.  lymphatic)  capillary 
wall  per  unit  volume  of  tissue.  In  the  special  case  y  =  0,  we  recover  the  linear  biphasic 
equations  describing  the  deformation  of  avascular  cartilage  like  materials. 

It  is  implicit  in  the  Equations  (1)  and  (2)  that,  in  general,  mechanical  loading 
not  only  strains  the  tissue,  but  also  pressurizes  both  solid  and  fluid  phases.  The 
pressurization  mechanism  can  be  understood  based  on  the  mechanical  behavior  of  the 
drained  interstitium  (elastic  solid  matrix).  In  contrast  to  the  solid  phase,  the  solid 
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matrix  is  compressible,  i.e. ,  it  reduces  its  volume  when  a  mechanical  loading  is  applied 
by  reducing  its  pore  space  volume.  In  an  ideal  case,  where  there  is  no  saturating  fluid 
or  the  saturating  fluid  can  move  frictionless  within  the  pore  system  and  drain  freely  to 
the  vascular  compartment,  the  pore  system  would  shrink  instantaneously  in  response 
to  the  applied  loading.  However,  as  the  interstitial  fluid  face  resistance  to  percolate  and 
drain,  it  resists  to  the  pore  system  shrinking,  pressurizing  and  being  pressurized  by  the 
solid  phase.  Relaxation  takes  place,  i.e.,  pressure  drops  gradually,  as  the  interstitial 
fluid  percolates  or  drains  in  response  to  the  pressurization.  During  the  relaxation 
the  tissue  approaches  to  the  solid  matrix  static  equilibrium.  At  static  equilibrium, 
p  =  0,  and  the  mechanical  behavior  is  governed  by  the  solid  matrix  Lame  parameters. 
It  follows  from  the  assumption  that  both  solid  and  fluid  phases  are  incompressible 
that  infinitesimal  dilatation  can  occur  only  when  the  corresponding  volume  of  fluid 
percolates  or  drains  to  vascular  or  lymphatic  systems. 

2.1.1.  One  dimensional  analytical  solution  Solutions  of  equations  (1)  and  (2)  exhibit 
various  stress/strain  relaxation  behaviors.  Insight  into  these  behaviors  and  their 
connection  to  different  physical  constants  that  appear  in  the  governing  equations  may 
be  developed  by  examining  a  simple  analytical  solution.  This  analytical  solution  also 
serves  as  a  check  on  the  validity  of  our  finite  element  (FEM)  numerical  implementation. 

For  the  purposes  of  developing  an  analytical  benchmark  solution,  we  consider  a 
2D  plane  strain  scenario  of  a  rectangular  homogeneous  tissue  sample  in  an  unconfined 
compression  test.  The  configuration  is  shown  schematically  in  Figure  1.  The  sample 
has  dimensions  of  L  x  h  and  fluid  can  flow  freely  across  the  lateral  boundaries  (p  =  0), 
which  are  also  traction  free.  The  fluid  cannot  flow  across  the  top  and  bottom 
boundaries,  which  are  also  shear  stress  free  (slip  boundary  conditions).  uy  =  0  at 
the  bottom  while  a  displacement  step  (or  ramp)  function  is  applied  at  the  top. 

For  a  step  function  compression  of  magnitude  uq,  the  pressure  field  in  the  sample 
is  given  by  (equation  (B.5)  in  Appendix  B) 

p{x,  t)  =  yy  exp(— x(A  +  2p)  t) 

00  3 

x  V (  — )  sin(a„;r)  exp(— cL?k(A  +  2p)  t)  (3) 

“  an 

n—  1 

77  7T 

Here,  an  =  — ,  0n  =  1  -  (-l)n. 

Note  that  the  pressure  relaxation  associated  with  microvascular  filtration, 
exp(— y(A  +  2p)t),  is  uniform  over  the  entire  sample.  This  is  the  result  that  a 
uniform  vascular  distribution  (as  assumed  in  this  example)  drains  all  parts  of  the 
tissue  at  the  same  rate.  The  pressure  relaxation  associated  with  fluid  percolation 
(the  exp(— a^«(A  +  2p)  t)  factor),  on  the  other  hand,  is  nonuniform  over  the  sample. 
The  exponential  factor  is  different  for  each  n  in  the  sum,  and  each  factor  multiplies 
a  spatial  “mode  shape”  sin(a„a;).  Thus  each  spatial  mode  decays  at  a  different  rate. 
While  the  mode  shapes  in  this  example  are  relatively  simple,  in  general  they  depend 
on  the  sample  geometry  and  boundary  conditions. 

The  horizontal  normal  strain,  exx,  behaves  very  similarly  to  the  pressure  in  this 
example.  Equation  (A. 3)  shows 

(A  T  2p)exx  —  p  Ae^y 

=  -  Ay  +p(x,t) 


(4) 

(5) 
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Thus  in  this  simple  example,  the  spatiotemporal  behavior  of  the  interstitial  pressure 
p  is  reflected  directly  in  the  lateral  strain  distribution. 

The  above  expressions  are  valid  for  a  step  compression  applied  instaneously  with 
time.  The  principal  of  linear  superposition  allows  us  to  integrate  in  time  equation  (5) 
to  obtain  the  response  for  a  gradual  compression  of  the  sample.  In  particular,  (see 
Appendix  B  for  details)  for  an  overall  compression  uo  applied  linearly  with  time  over 
an  interval  of  time  ti,  the  pressure  distribution  is  very  similar  in  form  to  (3);  viz. 


p(x,t) 


A 


n 


exp(— x(A  +  2 p)  t) 

oo 

x  An  sin(a„a;)  exp(— +  2p)  t) 

71=1 

/3»[exp ((q£k  +  y)(A  +  2p)t1)  -  1] 
an{a2nn  +  x)(A  +  2p)) 


for  t>t\.  (6) 

(7) 


The  lateral  strain  is  again  given  in  terms  of  the  pressure  by  equation  (5). 

Equations  (6)  and  (5)  indicate  that  the  strain  in  ^-direction  ( exx )  reaches  its 
maximum  right  after  the  mechanical  loading  has  been  applied  and  then  decreases 
with  time,  as  the  tissue  relaxes.  This  is  a  typical  pattern  in  unconfined  tests  of 
poroelastic  samples,  and  can  be  understood  by  considering  the  configuration  assumed 
by  the  sample  in  the  two  limits,  t  =  0+  and  t  —  oo.  Immediately  after  a  step 
mechanical  loading  has  been  applied  ( t  =  0+),  fluid  has  not  yet  had  a  chance  to  leave 
the  sample.  Thus,  at  t  =  0+,  we  expect  the  sample  to  behave  like  an  incompressible 
solid,  i.e.,  it  occupies  the  same  volume  (area)  as  it  occupied  previously.  On  the  other 
hand,  when  the  tissue  is  completely  relaxed  (<  =  oo),  the  pressure  relaxes  to  zero 
and  the  mechanical  behavior  is  governed  by  the  solid  matrix  Lame  parameters.  Since 
the  solid  matrix  is  compressible,  the  sample  should  now  occupy  a  smaller  volume 
(area).  As  uy  is  constant  with  time  in  this  experiment,  (due  to  the  constant  boundary 
conditions),  any  volume  reduction  must  be  reflected  in  a  shrinking  in  the  x  direction. 
Of  course,  during  the  transient  regime,  while  the  fluid  is  exuding  and/or  draining,  all 
the  configurations  assumed  by  the  sample  correspond  to  configurations  between  these 
two  extremes. 

As  mentioned  above,  in  both  equations  (3)  and  (6)  we  can  identify  in  two  different 
transient  phenomena,  percolation  and  vascular  drainage.  These  are  controlled  by 
«^k( A  +  2/i)  and  x(^  +  2 p),  respectively.  The  form  of  these  constants  indicates 
that,  beside  the  dependence  on  the  interstitial  hydraulic  conductivity  and  filtration 
coefficient,  the  greater  A  and  p  are  (especially  A  which  usually  is  much  larger  than  p), 
the  faster  the  tissue  relaxes.  The  ratio  of  k/ L 2  and  x  indicates  the  relative  importance 
of  fluid  flow  within  the  extravascular  compartment  versus  fluid  exchange  between 
compartments. 

To  illustrate  this,  we  plot  in  Figure  2  the  solution  for  Equation  (B.7),  at  three 
different  times,  for  L  =  10 an,  h  =  10 cm,  uq  =  1  cm  and  t\  =  0.3sec.  For  the 
blue  line,  the  poroelastic  parameters  are  chosen  such  that  vascular  drainage,  or  fluid 
exchange  between  compartments,  is  the  dominant  phenomenon.  (The  parameters  used 
are  listed  in  Table  I  and  correspond  to  the  inclusion  in  Experiment  1.)  We  observe 
that  the  pressure  is  homogeneous  along  almost  the  whole  sample  width.  The  pressure 
decreases  very  rapidly  close  to  the  lateral  surfaces,  which  indicates  that  only  a  small 
amount  of  fluid  crosses  the  boundary.  For  the  red  line,  we  increase  k  10,000  times  and 
decrease  x  1,000  times.  Now,  we  see  that  the  pressure  boundary  layer  rapidly  becomes 
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thicker  with  time,  indicating  that  percolation  is  now  dominant.  In  Figure  3  we  plot 
the  strain  field,  exx,  at  t  =  5.4sec,  for  both  cases.  In  the  Figure  3(a),  in  agreement 
with  the  pressure  field  plotted  in  Figure  2  (blue  solid  line),  we  see  that  the  strain  is 
approximately  constant  along  almost  the  whole  sample.  Again,  this  indicates  that 
the  vascular  drainage  is  dominant  and  percolation  is  negligible  outside  the  very  thin 
boundary  layers  located  at  the  sample  laterals.  It  means  that  exx  is  decreasing  with 
time  uniformly  with  x.  In  Figure  3(b)  we  see  that  the  strain  varies  with  x,  indicating 
that  percolation  is  now  important.  In  agreement  with  [3],  it  indicates  that  now  the 
shrinking  in  the  x  direction  diffuses  with  time  from  the  laterals  toward  the  center.  It 
is  interesting  to  note  that  at  t  =  0+  and  t  =  oo  the  sample  configurations  are  the  same 
in  both  cases,  the  difference  being  in  how  the  samples  go  from  one  limit  to  the  other. 

2.1.2.  Selection  of  material  properties  The  poroelastic  medium  is  defined  by  5 
physical  parameters:  <f>f  A,  p,  k  and  y.  We  work  in  plane  strain  state.  The  2- 
D  scenario  assumed  in  the  computational  simulations  is  that  of  a  circular  inclusion 
embedded  in  a  square  tissue  sample,  as  shown  in  Figures  4  and  5.  The  circular 
inclusion  is  intended  to  represent  a  malignant  breast  tumor,  with  stiffness,  microvessel 
density,  hydraulic  conductivity,  and  connective  tissue  density  all  elevated  relative  to 
the  background,  normal,  values. 

Normal  tissue  properties:  The  tissue  surrounding  the  inclusion  is  assumed  to  have 
“normal”  biomechanical  properties.  The  shear  modulus,  /i,  was  chosen  according  the 
values  reported  for  breast  tissue.  We  then  calculated  the  corresponding  value  for  A 
assuming  a  Poisson  ratio  of  0.49.  The  values  for  (f>f ,  n  and  \  were  chosen  according 
to  values  given  in  [6].  In  Experiment  3,  we  arbitrarily  increased  the  background 
interstitial  hydraulic  conductivity,  n ,  to  investigate  percolation  effects. 

Malignant  tissue  properties:  Regarding  the  inclusion,  in  order  to  reproduce  a  solid 
tumor  behavior,  we  assumed  an  augmented  capillary  filtration  coefficient  [6]  in  all  the 
experiments.  We  also  assumed  it  is  stiffer,  increasing  the  value  of  p1  and  assumed 
a  Poisson  ratio  of  0.47,  calculating  A  accordingly.  All  the  parameters  used  are 
summarized  in  Table  I. 

2.2.  Computational  experiments 

In  order  to  evaluate  the  predictions  of  the  mathematical  model  with  nontrivial 
geometries  and  boundary  conditions,  we  developed  a  finite  element  discretization 
of  equations  (1)  and  (2)  in  two  dimensions.  We  used  the  standard  Galerkin 
approximation  with  bilinear  shape  functions  for  both  the  pressure  and  displacement 
fields.  To  integrate  in  time  we  use  the  Backward  Euler  method,  assuming  all 
material  parameters  are  constant  with  time.  We  have  validated  our  implementation 
by  comparing  the  numerical  solution  to  the  analytical  solution  derived  in  the  previous 
section,  as  shown  in  the  Figure  2. 

We  now  use  this  finite  element  implementation  to  study  two  dimensional  problems 
that  model  hypothetical  clinical  imaging  exams.  In  the  computational  experiments 
presented  here,  we  tried  to  reproduce  hypothetical  configurations  for  clinical  breast 
exams,  and  investigate  the  strain  relaxation  within  the  sample. 


Interstitial  Fluid  Flow 


39 


2.2.1.  Experiment  1  Experiment  1  is  schematically  shown  in  Figure  4.  The  circular 
inclusion  has  1  cm  diameter  and  the  sample  has  dimensions  10cm  x  10cm.  The  fluid 
cannot  flow  across  the  boundaries,  mimicking  a  portion  of  tissue  completely  bounded 
by  skin.  Therefore,  the  interstitial  fluid  can  redistribute,  but  the  only  way  for  it  to 
leave  the  sample  is  by  vascular  drainage.  Such  an  idealized  boundary  condition  is 
valid  when  the  drainage  effects  are  much  larger  than  the  percolation  effects  and  the 
permeable  boundary  is  relatively  far  away  from  the  region  to  be  investigated,  as  in  the 
case  here.  The  tissue  is  fixed  at  the  bottom,  where  ux  =  uy  =  0.  The  lateral  surfaces 
are  traction  free.  At  the  top,  we  simulate  the  mechanical  loading  from  a  compressor 
of  5  cm.  of  width.  The  displacement  of  the  compressor  is  modeled  by  a  ramp  function 
such  that  the  prescribed  uy  goes  from  0  to  1  cm.  in  0.3  sec,  in  the  region  corresponding 
to  x  =  2.5 cm  to  x  =  7.5cm.  Below  the  compressor  we  prescribe  zero  shear  stress 
(tvx),  which  models  a  slip  boundary  condition. 

2.2.2.  Experiment  2  Experiment  2  is  schematically  shown  in  Figure  5.  The  circular 
inclusion  still  has  1cm  diameter,  while  the  model  still  has  dimensions  of  10  cm  x  10  cm. 
As  before,  the  fluid  cannot  flow  across  the  boundaries  and  the  tissue  is  fixed  at  the 
bottom.  However,  now  the  model  is  completely  confined  at  the  top,  where  uy  goes  from 
0  to  0.03cm  in  0.3sec,  and  is  partially  confined  at  the  lateral  surfaces,  i.e.,  ux  =  0  from 
y  =  2.0 cm  to  y  =  10cm.,  while  it  is  traction  free  from  y  =  0.0 cm.  to  y  =  2.0cm.  The 
goal  here  is  to  reproduce  a  situation  of  partial  breast  confinement,  with  the  recognition 
that  the  breast  cannot  be  completely  confined  in  the  clinic. 

2.2.3.  Experiment  3  Experiment  3  has  the  same  configuration  as  Experiment  1,  but 
the  interstitial  permeability,  k,  is  the  same  for  both  inclusion  and  surrounding  tissue 
and  is  6.4  x  10~13m2/(Pa.sec). 

2.2.4-  Experiment  4  Experiment  4  has  the  same  configuration  as  Experiment  1, 
however  the  fluid  can  flow  freely  across  the  lateral  surfaces  of  the  sample,  where  p=0. 

3.  Results 

In  this  section  we  show  results  corresponding  to  the  region  delimited  by  the  dotted 
line  in  Figures  4  and  5.  It  has  dimensions  of  4  cm  x  4  cm  and  is  contained  between 
x  =  3.0  cm  to  x  =  7.0  cm  and  y  =  5.5 cm  to  y  —  9.5cm.  We  emphasize  that  we  have 
solved  the  problem  in  the  entire  domain,  but  are  showing  the  results  only  in  this  region 
of  interest,  in  order  to  investigate  the  behavior  of  the  inclusion  and  its  surroundings  in 
detail.  This  is  intended  to  be  representative  of  ultrasound  imaging  where  the  physical 
boundaries  are  typically  distinct  from  the  image  boundaries. 

3.1.  Experiment  1 

In  this  experiment,  the  fluid  exchange  between  interstitial  and  microvascular 
compartments  is  the  dominant  phenomenon.  Due  to  the  difference  between  the 
filtration  coefficient  inside  and  outside  the  inclusion  (it  is  30  times  larger  inside), 
a  transient  analysis  of  the  problem  can  be  outlined  considering  two  different  time 
scales:  the  inclusion’s  relatively  short  relaxation  time  and  the  surrounding  tissue’s 
large  relaxation  time. 


Interstitial  Fluid  Flow 


40 


Right  after  the  mechanical  loading  has  been  applied  and  before  significant 
fluid  drainage  has  occurred,  the  tissue  is  pressurized  and  the  sample  approximately 
behaves  like  an  incompressible  elastic  solid,  with  the  same  shear  modulus  (/z)  as 
the  corresponding  solid  matrix  and  a  Young’s  modulus  equals  to  3fi.  The  pressure 
held  at  t  =  0.3 sec  is  shown  in  Figure  6(a).  We  can  see  the  stress  concentrations  at 
the  transducer  edges  radiating  in  the  upper  left  and  right  corners  of  the  figure.  At 
the  center,  we  can  distinguish  the  inclusion  and  four  lobes  resulting  from  the  stress 
concentration  at  the  inclusion. 

Gradually,  as  the  fluid  drains  from  the  interstitium  to  the  microvasculature, 
the  tissue  relaxes.  In  Figure  6(b)  we  plot  the  pressure  held  at  t  =  5.4 sec. 
Comparing  Figures  6(a)  and  6(b)  shows  that  the  inclusion  relaxes  much  faster 
than  the  surrounding  tissue.  This  may  be  attributed  to  the  higher  value  of  the 
microvascular  filtration  coefficient  in  the  inclusion.  In  Figure  6(c)  we  plot  the 
difference  in  the  x-normal  strain  held  (exx)  between  t  =  0.3sec  and  t  =  5.4sec,  i.e., 
exx(t-5Asec)  —  exx(t=o.3sec)-  In  Figure  6(d)  we  plot  the  difference  in  the  y-normal  strain 
held,  i.e.,  eyy(t=5.4Sec)  -  eyy(t=o.3sec )  and  in  Figure  6(e)  we  plot  the  corresponding 
dilatation  difference,  i.e.,  (ex'2:(z— 5.4*ec)  5.4sec) )  (Gx(z— o.3sec)“b^yy(t=o.3sec))-  ^Ye 

observe  that  the  dilatation  in  the  inclusion  is  negative,  indicating  that  it  is  shrinking 
as  it  relaxes.  On  the  other  hand,  the  volume  of  the  surrounding  tissue  remains  almost 
unchanged  (dilatation  «  0).  We  also  observe  four  lobes  in  Figures  6(c)  and  6(d) 
resulting  from  the  stress  concentration  redistribution  around  the  inclusion,  which  has 
occurred  during  the  relaxation.  We  can  see  in  Figure  6(e)  that  the  huid  drainage  at 
the  lobes  is  small,  since  the  dilatation  in  this  region  is  almost  zero.  In  summary,  at 
t  =  5.4sec,  the  mechanical  behavior  of  the  inclusion  is  approximately  governed  by  the 
corresponding  solid  matrix  mechanical  properties;  that  is,  it  has  relaxed,  while  the 
surrounding  tissue  still  behaves  like  an  incompressible  material. 

The  inclusion  takes  about  10  secs  to  relax  almost  completely.  In  the  Figure 
6(f)  we  plot  the  pressure  Held  at  t  =  10.2  sec.  We  observe  that  the  pressure  inside 
the  inclusion  is  close  to  the  equilibrium  pressure  ( p  «  0).  In  Figure  6(g)  we  plot 
Gs(t=i5.0sec)  -  Gx(t=io.2Sec),  in  Figure  6(h)  we  plot  eyy(t  =  15.0sec)  €yy(t=10.2sec)  Slid 

in  Figure  6(i)  we  plot  (eXa;(t=15.0 sec)  +  eyy(t=15.0sec))  ~  (exx(t=  10.2sec)  +  £yy(t=10.2sec))- 
Now,  both  the  inclusion  and  the  surrounding  tissue  relax  by  similar  amounts.  Due 
to  the  applied  displacement,  we  see  that  the  strain  occurs  predominantly  in  the  x 
direction  and  both  inclusion  and  surrounding  tissue  shrink  at  approximately  the  same 
rate.  We  also  observe  that  the  interstitial  fluid  drains  faster  (or  percolates)  in  a  thin 
region  around  the  inclusion,  due  to  the  stress  concentration.  The  surrounding  tissue 
takes  about  300sec  to  relax  almost  completely.  As  discussed  before,  at  the  steady 
state,  where  both  inclusion  and  surrounding  tissue  are  relaxed,  the  sample  assumes 
the  configuration  where  the  mechanical  behavior  of  both  inclusion  and  surrounding 
tissue  are  governed  by  the  respective  solid  matrix  Lame  parameters. 

3.2.  Experiment  2 

Here,  as  before,  the  fluid  exchange  between  interstitial  and  microvascular 
compartments  is  the  dominant  phenomenon.  A  transient  analysis  of  the  problem 
can  again  be  outlined  by  considering  two  different  time  scales.  The  pressure  field  at 
t  =  0.3 sec  is  shown  in  Figure  7(a).  We  see  that  the  pressure  magnitude  is  similar  to  the 
previous  case  despite  the  much  smaller  displacement  prescribed  at  the  top  boundary. 
This  is  obtained  by  confining  a  portion  of  the  lateral  surface.  We  also  observe  that 
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in  contrast  to  the  previous  case,  there  is  no  stress  concentration  in  the  upper  left  and 
right  corners  and  the  pressure  field  is  almost  uniform. 

As  the  interstitial  fluid  drains,  the  inclusion  relaxes.  The  pressure  field  at 
t  =  6.3 sec  is  shown  in  Figure  7(b).  We  can  see  that  the  pressure  inside  the  inclusion  is 
about  85%  smaller  at  t  =  6.3 sec  than  at  t  =  0.3 sec,  while  it  remains  close  to  the  peak 
in  the  surrounding  tissue.  Once  again  this  is  due  to  the  higher  vascular  filtration  in  the 
inclusion.  In  the  Figure  7(c)  we  plot  exx(t=6.3sec)  -eXx(t=o.3«ec)>  in  Figure  7(d)  we  plot 
e yy(t—6.3sec )  ^yy(t=0.3sec)  ^rid  in  in  Figure  /(e)  we  plot  (exx(i=6.3sec)  T  ^yy{t— 6 .3sec)) 

(eXx(t=o.3sec)  +  eyy(t=o.3sec))-  As  in  the  previous  case,  we  observe  that  the  dilatation  in 
the  inclusion  is  negative,  indicating  that  it  is  shrinking  as  it  relaxes,  and  the  volume 
of  the  surrounding  tissue  remains  almost  unchanged  (dilatation  «  0).we  also  see  four 
lobes  resulting  from  the  stress  concentration  redistribution  around  the  inclusion. 

The  inclusion  takes  about  10  sec  to  relax  almost  completely.  The  pressure 
field  for  t  =  13.5  sec  is  shown  in  Figure  7(f).  The  pressure  inside  the  inclusion  is 
now  negative,  indicating  that  fluid  is  draining  back  to  the  interstitial  compartment 
from  the  microvasculature,  as  the  surrounding  tissue  relaxes.  In  Figures  7(g), 
7(h)  and  7(i)  we  plot  exx(£=3l.5sec)  ^xx(i=13 .5sec),  Cyy(t= 31.5sec)  ^yy(t=13.bsec )  &rid 
iexx(t=3i.5sec)  +  eyy(t—3i.5sec))  ~  (e*x(t=i3.5sec)  +  eyy(t-i3.5sec)) i  respectively.  In  contrast 
to  the  previous  case,  due  to  the  boundary  conditions,  the  strain  predominantly  occurs 
in  the  y  direction.  Further,  in  agreement  with  the  pressure  field  shown  in  the  Figure 
7(f),  we  observe  that  during  this  period  the  inclusion  swells  while  the  surrounding 
tissue  shrinks.  In  fact,  the  inclusion  attains  its  smallest  volume  at  approximately 
t  =  13.5  and  then  swells  till  the  sample  reaches  the  steady  state  configuration. 

3.3.  Experiment  3 

In  this  experiment  we  investigate  the  impact  of  increasing  the  interstitial  permeability, 
on  the  spatio-temporal  strain  pattern.  The  permeability  here  is  100  times  higher 
than  in  the  previous  cases,  and  is  the  same  for  both  background  and  inclusion. 
Here  the  percolation  due  to  the  interstitial  pressure  gradient  at  the  periphery  of  the 
inclusion  is  significant.  The  pressure  field  at  t  =  0.3sec  is  shown  in  Figure  8(a). 
We  observe  that  it  is  similar  to  the  one  shown  in  the  Figure  6(a)  albeit  it  is  more 
diffuse.  In  the  Figure  8(b)  we  plot  the  pressure  field  at  t  =  5.4 sec.  Differently  from 
before,  we  can  distinguish  an  elliptic  region  at  the  inclusion’s  periphery,  where  the 
pressure  decreases  slower  then  the  central  region,  indicating  that  fluid  flowed  from 
the  surrounding  tissue  toward  the  inclusion,  partially  replacing  the  drained  fluid.  In 
Figures  8(c),  8(d)  and  8(e)  we  plot  Cxx(£=5_4sec)  £xx(f=0.3secp  ^ yy(t=5Asec )  ^yy(t=0.3sec) 

and  (eXa,(t=5.4sec)  ~^~^yy(t=5.4sec) )  (^xx(t=o.3sec) ~I^yy(t=o.3sec)) >  respectively.  In  contrast 
to  Experiment  1,  we  can  see  four  lobes  in  the  dilatation  field,  confirming  that  some 
fluid  has  percolated  around  the  inclusion  in  response  to  the  stress  concentration. 

In  the  Figure  8(f)  we  plot  the  pressure  field  for  t  =  10.2sec.  In  Figure  8(g), 

8(h)  and  8(i)  we  plot  Cxx(f=15.0sec)  ^xx(i=U0.2sec)>  ^yy(t=15.0sec)  ^ yy(t=10.2sec )  &nd 
(^xx(£=15.0sec)  T  ^yy(t=15.0sec))  (exx(t=io.2sec)  +  ej/y(t=io.2sec)))  respectively.  We  see 
that  the  thin  region  around  the  inclusion  identified  in  Figure  6(i)  is  now  (see  Figure 
8(i))  much  thicker.  The  sample  takes  about  300sec  to  relax  almost  completely.  At  the 
steady  state,  it  assumes  the  same  configuration  as  in  Experiment  1. 
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3-4-  Experiment  4 

This  experiment  is  identical  to  Experiment  1,  except  now  the  lateral  surfaces  are 
permeable  (p  =  0).  In  the  Figure  9(a),  (b)  and  (c)  we  plot  the  pressure  field  at 
t  =  0.3,  5.4  and  10.2sec,  respectively.  We  plot  the  solution  in  the  entire  domain 
because  the  only  change,  when  compared  with  experiment  1,  is  the  appearance  of 
thin  boundary  layers  along  the  lateral  surfaces.  Due  to  the  small  interstitial  hydraulic 
conductivity  the  effect  of  the  permeable  boundary  condition  is  confined  to  this  thin 
layer. 

4.  Discussion 

In  Experiment  1,  we  observe  that  in  the  first  5  seconds  the  inclusion  relaxes  much 
more  than  the  surrounding  tissue.  During  this  time,  its  volume  reduces  by  about 
4, 000  microstrains.  For  an  inclusion  of  1cm.  diameter  this  implies  displacements  of 
the  order  of  40  microns.  On  the  other  hand,  after  the  inclusion  stops  relaxing,  the 
surrounding  tissue  relaxation  becomes  the  predominant  phenomenon.  Since  drainage 
is  the  main  relaxation  mechanism  for  both  inclusion  and  surrounding,  their  relaxation 
time  scales  are  directly  proportional  to  the  respective  filtration  coefficient  values  and, 
as  seen  in  Equations  (B.5)-(B.8),  can  be  roughly  estimated  by  r  =  In  terms 

of  the  shrinking  magnitude  resulting  from  the  relaxation,  it  is  directly  proportional 
to  the  applied  displacement  at  the  boundary.  Further,  similarly  to  what  is  discussed 
in  Section  2.2,  it  is  also  directly  proportional  to  the  Poisson’s  ratio  difference  between 
incompressible  material  (0.5)  and  the  solid  matrix,  i.e. ,  the  smaller  the  solid  matrix 
Poisson’s  ratio,  the  larger  the  shrinking. 

As  in  the  previous  experiment,  we  see  in  Experiment  2  that  in  the  first  6  seconds 
the  inclusion  shrinks  about  4,  000  microstrains,  which,  again,  for  an  inclusion  of  1  cm 
diameter  implies  displacements  of  the  order  of  40  microns. The  applied  displacement, 
the  unconfined  to  confined  lateral  area  ratio  and  the  inclusion  solid  matrix  bulk 
modulus  determine  the  shrinking  magnitude.  In  the  limit  case,  where  the  sample 
is  completely  confined,  practically  all  the  sample’s  volume  reduction  resulting  from 
the  applied  displacement  must  be  reflected  in  the  inclusion’s  volume  reduction.  We 
observed  also  that  after  the  inclusion  relaxes  it  experiences  a  gradual  swelling,  while 
the  surrounding  relaxation  takes  place.  It  can  be  understood  by  recognizing  that  the 
surrounding  tissue  shrinks  as  it  relaxes,  because  its  solid  matrix  is  compressible.  The 
shrinking  is  partially  balanced  by  the  inclusion’s  swelling.  Again,  in  the  limit  case, 
where  the  sample  is  completely  confined,  all  the  shrinking  must  be  balanced  by  the 
inclusion’s  swelling. 

In  Experiment  3,  we  can  observe  that  raising  the  interstitial  permeability  by 
a  factor  of  100  had  a  very  small  effect  on  the  spatio-temporal  patterns  of  strain, 
especially  in  the  short  term.  In  a  comparison  with  Experiment  1,  we  observe  that 
the  pressure  field  is  no  longer  homogeneous  inside  the  inclusion,  decreasing  from  the 
periphery  toward  the  center,  indicating  that  part  of  the  drained  fluid  is  being  replaced 
(at  the  periphery)  by  fluid  flowing  from  the  surrounding.  For  the  same  reason,  part 
of  the  surrounding  tissue  close  to  the  inclusion  experiences  now  a  larger  dilatation.  In 
the  middle  and  long  terms  the  impact  is  larger,  especially  around  the  inclusion.  There 
we  could  see  the  dilatation  partially  “diffusing”  from  the  center  of  the  sample  toward 
the  boundaries.  Here,  an  analogy  can  be  made  with  what  is  discussed  in  Section 
2.2,  where  now  the  inclusion’s  limit  plays  the  role  of  the  open  lateral  boundaries  in 
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Figure  2.  In  a  case  where  the  interstitial  permeability  is  as  important  as  (or  more 
important  than)  the  filtration  coefficient,  a  transient  analysis  of  the  problem  could 
no  longer  be  outlined  considering  two  different  time  scales  -  the  inclusion’s  relatively 
short  relaxation  time  and  the  surrounding  tissue’s  large  relaxation  time  -  despite  the 
difference  between  the  filtration  coefficient  inside  and  outside  the  inclusion.  The  whole 
sample  would  relax  approximately  with  the  same  rate  due  to  the  fluid  replacement 
mechanism  discussed  above. 

The  effect  of  letting  the  fluid  freely  flow  across  the  lateral  boundaries  in 
Experiment  4  is  imperceptible  inside  the  region  of  interest  investigated  in  the  previous 
experiments.  In  fact,  its  effect  is  only  felt  in  a  very  thin  region  close  to  the  sample’s 
laterals,  similarly  to  what  is  discussed  in  Section  2.2.  This  is,  of  course,  consistent 
with  the  notion  that  the  fluid  exchange  between  compartments  is  much  more  important 
than  the  fluid  flow  within  the  extravascular  compartment.  However,  in  cases  where 
the  interstitial  permeability  is  significant,  permeable  boundaries  become  important. 

The  results  suggest  that  it  may  be  possible  to  image  the  interstitial  fluid  motion 
in  tissues  by  measuring  the  corresponding  strain  rate.  A  sequence  of  images  acquired 
from  ultrasound  or  other  scanners  could  be  processed,  as  they  are  in  elasticity  imaging, 
to  track  the  spatio-temporal  patterns  of  elastic  strain.  In  addition,  the  strain  pattern 
could  then  be  used  to  solve  for  the  spatial  distribution  of  the  poroelastic  parameters, 
in  particular,  the  shear  modulus  /i  and  the  microvascular  filtration  coefficient  \- 

It  is  interesting  to  consider  the  ultrasound  measurability  of  the  transient  strains 
predicted  here.  In  experiment  1,  we  noted  a  volume  change  in  the  inclusion  of  about 
0.4%  after  5  seconds,  in  an  overall  compression  of  10%.  Such  a  relaxation  would 
certainly  be  measurable  by  ultrasound,  though  tracking  a  compression  over  a  full  10% 
strain  might  present  technical  difficulties.  On  the  other  hand,  in  experiment  2,  with 
confined  compression  of  0.3%  we  noted  the  same  inclusion  volume  change  of  0.4% 
over  about  6  seconds.  The  volume  change  is  roughly  isotropic  in  the  plane,  so  about 
half  that  of  that,  or  about  0.2%,  would  take  place  in  the  high  resolution  (ultrasound 
propagation)  direction.  In  practice,  it’s  likely  that  the  plane  strain  assumption  would 
be  violated  here,  and  the  volume  change  might  be  expected  to  be  isotropic  in  the 
volume.  In  that  case,  only  about  1/3  of  the  total  volume  change,  or  about  0.15% 
strain,  would  be  reflected  in  the  high  resolution  direction.  This  magnitude  of  strain, 
over  6  seconds,  would  certainly  be  a  measurable  effect. 

Whether  the  magnitudes  of  the  effects  predicted  here  might  be  seen  in  clinical 
practice  depends  on  the  validity  of  several  assumptions.  First,  these  magnitudes 
depend  on  the  validity  of  the  parameters  chosen  for  the  model.  These  are  selected 
as  described  above,  and  may  be  assumed  to  have  large  variability  in  practice.  Our 
model  assumes  a  linear  response  to  the  applied  pressure.  This  implies,  among 
other  limiting  situations,  that  neither  the  porosity  nor  permeability  will  change  with 
applied  pressure.  It  is  likely  that  this  second  order  effect  will  indeed  remain  small 
provided  the  applied  pressure  in  the  tissue  stays  well  below  the  collapse  pressure  of 
the  microvasculature.  The  model  neglects  all  forms  of  transient  solid  response  that  is 
not  directly  related  to  fluid  flow.  The  magnitude  of  these  possible  effects  is  impossible 
to  estimate  at  this  time.  Finally,  the  plane  strain  assumption  is  taken  as  analytically 
convenient  for  our  purposes,  but  is  not  supposed  that  it  is  quantitatively  accurate. 
Quantitative  discrepancies  between  2 D  and  3 D  predictions  of  as  much  as  50%  might 
reasonably  be  expected,  though  order  of  magnitude  changes  would  not  be  expected. 
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5.  Conclusions 

A  poroelastic  model  that  includes  the  effects  of  fluid  flow  and  the  possibility  of 
exchange  between  fluid  compartments  was  used  in  conjunction  with  finite  element 
modeling  to  predict  the  effects  of  fluid  flow  on  the  spatio-temporal  patterns  of  soft- 
tissue  elastic  strain  under  a  variety  of  physiological  conditions. 

The  analytical  solution  for  the  problem  consisting  of  a  tissue  sample  in  an 
unconfined  test  show  two  different  transient  phenomena,  percolation  and  drainage, 
controlled  by  k  and  respectively.  In  soft  tissues  drainage  is  the  dominant 
phenomenon. 

Numerical  simulation  results  suggest  that  it  may  be  possible  to  image  the 
interstitial  fluid  motion  in  tissues  by  measuring  the  corresponding  strain  rate. 
Further,  they  show  that  the  abnormal  tumor  microvasculature  may  increase  the  strain 
relaxation  rate. 

In  unconfined  tests  the  total  dilatation  resulting  from  the  tissue  relaxation  is 
controlled  by  the  Poisson’s  ratio  difference  between  incompressible  material  (0.5)  and 
the  solid  matrix,  while  in  partially  confined  tests  it  is  controlled  by  the  unconfined 
to  confined  lateral  area  ratio.  We  didn’t  perturb  significantly  the  spatio-temporal 
patterns  of  strain  by  increasing  the  interstitial  permeability  by  100  times,  specially 
while  the  inclusion  relaxes.  The  effect  of  letting  the  fluid  freely  flow  across  the  lateral 
boundaries  is  confined  to  a  very  thin  region  close  to  the  sample’s  laterals. 
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Appendix  A.  Derivation  of  field  equations 

In  this  section  we  present  a  rederivation  of  the  field  equations  from  [5].  Figure  10 
represents  schematically  a  portion  of  soft  tissue  interstitium.  Its  boundary  and  domain 
are  denoted  by  T  and  ft,  respectively.  The  interstitial  boundary,  T  is  comprised 
of  three  parts:  the  outer  boundary  F0,  the  interface  with  the  hemal  capillaries 
Fc,  and  the  interface  with  the  lymphatic  capillaries  Fl-  We  regard  the  interstitial 
compartment  as  a  linear  biphasic  solid-fluid  mixture,  where  both  fluid  and  solid  phases 
can  move  independently  of  each  other.  The  two  phases  are  treated  as  intrinsically 
incompressible.  Thus  the  total  stress  in  interstitium,  eU,  is  given  by: 

^  =  ~pl  +  XV -ul +  2iiVSymu  (A.l) 

Here  A  and  p  are  the  solid  matrix  Lame  parameters,  p  is  the  interstitial  fluid  pressure, 
and  u  is  the  displacement  of  the  solid  phase.  The  infinitesimal  strain  tensor  is  given 
by  e  =  VSymu  =  ^(Vu+(Vu)T).  Note  that,  in  contrast  to  the  solid  phase,  the  solid 
matrix  is  compressible. 

For  frequencies  and  rates  of  strain  which  occur  physiologically,  and  under  many 
clinical  applications,  the  effects  of  inertia  in  tissue  dynamics  may  be  neglected.  Under 
this  assumption  and  in  the  absence  of  body  forces,  the  resulting  momentum  equations 
for  the  fluid  phase  and  mixture  may  be  written  as  ([3]  and  [7]): 

V  ■  cr^  +  rF1w  =  0 

V  •  <r‘  =0 


(A.2) 

(A.3) 
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Here,  cr?  =  — pi  is  proportional  to  the  stress  in  the  fluid  phase,  k  is  the  interstitial 
permeability,  and  w  is  volume  average  relative  fluid  displacement  (relative  to  the 
solid  matrix).  Equation  (A. 2)  is  a  statement  of  Darcy’s  law  for  the  fluid  flow  through 
the  interstitium,  while  equation  (A. 3)  is  a  statement  of  conservation  of  total  linear 
momentum. 

For  a  biphasic  material  with  incompressible  phases,  the  equations  (A.l),  (A. 2) 
and  (A. 3)  are  augmented  by  the  incompressibility  constraint  for  the  fluid  phase: 

V  •  ii  +  V  •  w  =  0  (A. 4) 

To  bring  out  explicitly  the  effects  of  the  vascularization,  we  shall  average  (A. 4)  over  an 
elementary  “averaging  volume”,  that  is,  a  small  volume  (typically  0(lmm3)),  which 
is  large  enough  to  contain  a  sufficiently  large  number  of  microvessels  that  the  averages 
below  become  sensibly  stationary. 

To  that  end,  we  integrate  (A. 4)  over  Cl,  our  averaging  volume: 

I  V  •  u dV  +  I  V-wdV  =  0.  (A. 5) 

Jn  Jn 

We  apply  the  divergence  theorem  to  the  second  term  in  equation  (A. 5)  to  obtain: 


I  \/-udV  +  I  n  ■  w  dS  +  f  n  ■  w  dS  +  f  n  ■  w  dS  =  0.  (A. 6) 

Jn  Jr0  Jvc 

Here,  we  recall  that  divided  the  total  boundary  of  the  interstitial  space  into  three 
parts,  the  external  boundary  To,  the  (hemal)  capillary  surface  Tc,  and  the  lymphatic 
capillary  surface,  T^. 

The  last  two  terms  in  equation  (A. 6)  represent  fluid  fluxes  into  the  microvessels. 
As  such,  they  can  be  related  to  the  pressure  difference  across  the  microvessel  wall 
through  Starling’s  law  applied  to  both  classes  of  microvessels: 


Jc  =  Lp(pv  —  p)  (A. 7) 

Jl  =  Lpl(pl  —  p)  (A. 8) 

Here  J  is  the  volume  fluid  flux  out  of  the  vessel,  Lp  represents  the  hydraulic 
conductivity  of  the  microvessel  wall,  and  the  subscripts  c  and  L  denote  hemal 
capillaries  and  lympthatic  capillaries,  respectively.  pv  and  pl  thus  represent  the  hemal 
capillary  and  lymphatic  capillary  pressures,  respectively. 

Using  equations  (A. 7)  and  (A. 8)  allows  us  to  make  the  following  approximations: 


n  ■  w  dS  «  LpSv  ( p  —  pv) 


(A.9) 


/  n  ■  wdS  «  LpLSl  {p  —  Pl)  (A. 10) 

JrL 

Here,  Sv  and  Sl  are  the  total  surface  areas  of  the  two  classes  of  microvessels  within 
the  averaging  volume,  Cl. 

Further,  ensemble  averaging  the  first  and  second  terms  of  (A. 6)  over  all 
representative  realizations  of  the  microvasculature  yields: 


/  V  •  iidV  s 
Jn 

s  DV-  <  u> 

(A.ll) 

/  n  ■  wdS  p 

s  /  n-  <  w  >  dS  = 

/  V-  <  w  >  dV  =  DV-  <(&.a>2) 

1  r0 

J  r0 

Jn 
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Angle  brackets  in  (A. 11)  and  (A. 12)  represent  ensemble  averages.  We  now  drop  the 
angle  brackets  with  the  understanding  that  all  symbols  represent  quantities  averaged 
over  the  microstructure.  V  is  the  volume  of  Q. 

We  now  use  (A.9-A.12)  in  (A. 5)  to  write: 

V  •  u  +  V  •  w  +  xv  {p  -  Pv)  +  Xl  (p  -  Pl ) 
pFv 

x,  =  — 

LplSl 
Xl  =  — y~ 

Equation  (A.  13)  is  identical  to  Equation  4  in  [5]  with  appropriate  reinterpretation 
of  average  fluid  and  solid  displacements.  At  normal  physiological  conditions,  where 
Pv  >  p  >  Pl,  Xv{pv  ~  p)  and  Xl(p  —  Pl)  represent  the  transcapillary  flow  and  the 
lymphatic  drainage,  respectively. 

Transcapillary  exchange  and  lymphatic  drainage  can  be  expected  to  be  taking 
place  continuously  in  living  tissues.  Under  these  “steady  state  operating  conditions,” 
p  and  w  will  in  general  be  nonzero.  Therefore,  we  let 


u 

=  uss  +  ii 

(A. 16) 

■w 

=  wss  +  w 

(A.17) 

p 

=  pss+p 

(A.18) 

Pv 

=  PvSS  +Pv 

(A.19) 

PL 

=  Pl  +Pl 

(A. 20) 

where,  uss  is  the  (temporally)  constant  steady  state  part  of  it,  and  u  is  its  (temporally) 
fluctuating  part. 

Then  equations  (A.l),  (A. 2),  (A. 3)  and  (A.  13)  remain  valid  with  u  replaced  by 
u,  p  replaced  by  p,  etc.  In  particular,  in  the  special  case  pv  =  Pl  =  0,  equation  (A. 13) 
becomes 

V  •  u  +  V  •  w  +  xp  =  0.  (A. 21) 

Here,  X  —  Xv  F  XL  is  the  total  microvascular  filtration  coefficient. 

Substituting  equation  (A. 21)  into  (A. 2)  and  dropping  the  tildes  yields  equation 
(1).  Substituting  equation  (A.l)  into  (A. 3),  using  equations  (A.16-A.20)  and  dropping 
the  tildes  gives  equation  (2). 


=  0  (A. 13) 

(A.14) 

(A.15) 


Appendix  B.  One  dimensional  analytical  solution 


Under  the  assumptions  outlined  in  section  2.1.1,  we  find  that  the  pressure  and 
lateral  displacement  are  both  independent  of  the  vertical  coordinate,  i.e.  p  =  p(x,  t) 
and  ux  =  ux{x,t).  Further,  the  lateral  normal  stress  vanishes,  i.e.  olxx  =  0. 
Finally,  for  a  step  function  applied  displacement,  the  vertical  displacement  is  simply 
uy  =  uv(y,t)  =  —  uo-ff(t)(f),  where  H(t)  is  the  step  function.  Thus  the  problem 
reduces  to  a  ID  problem  and  Equations  (??)  and  (??)  reduce  to: 


dux 
dx 

I-<A 


d2p 
dx2  ^ 

2p) 


XP 

d2ux 

dx2 


=  0 


=  0 


(B.l) 

(B.2) 
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Integrating  Equation  (B.2)  over  x  and  using  Equation  (??)  (note  that  '-^Pr  =  -^a), 


dy 


gives: 


/  \  i  o  \  d‘Ux  _  ,  M0 

p  -  (\  +  2p)——  —  —  X— 

ax  h 


(B.3) 


Taking  the  time  derivative  of  Equation  (B.3)  and  using  it  in  Equation  (B.l),  yields 


p~kU+kp  =  0 


(B-4) 


Here  we  have  introduced  the  symbols  k  =  k( A  +  2/z)  and  %  =  y(A  +  2 p). 

Equation  (B.4)  can  then  be  solved  using  separation  of  variables  to  obtain: 

p(x,t)  =  V'C  — )  sin(ana;)e_7”t  (B.5) 

Lh  (  a„ 

n— 1 

where  an  =  — ,  =  1  —  (—1)"  and  7ra  =  Using  Equation  (B.3)  and  the 

1j 

fact  that  ux  =  0  at  x  =  j,  for  all  time  gives: 


u  (x  t)  = _ ^U°L _ (x  -  —  ) 

x[,)  h(\  +  2p)[X  2 

4pu0 


-  Lh( A+V)  S  i*008'0"1*  -  c»(T))e"7"‘  <B'6) 

To  extend  the  solution  to  a  case  where  a  ramp  function  is  applied,  we  take 
advantage  of  the  fact  that  a  ramp  is  the  integral  of  the  step  function.  Therefore 
we  integrate  Equations  (B.5)  and  (B.6)  from  t  =  0  to  t  =  t\.  This  yields: 


P(x,t)  =  < 


and 


ux(x,t)  = 


4 an  ,  1  -  e  lnt 

—  >  ( —  )sm(avr) - 

Lh  ^  an  'in 

n—1 

4 .  ,  e7»(ti-t)  _  e-7»t 

—  >  ( — )  sm(a„x) - 

Lh  z—  an  7n 


Aa(:r  -  —  )t  -  Aln  ,t  <  t\ 

£ 

Xa(x  —  —)t  —  A2„  ,t>ti 


,t<h 


,  t  >  ti 


(B.7) 


where 


4/zd  ^  /?„,  nrr  1-e7"4 

Ai„  =  -7—  >  —  (cos(a„:r)  -  cos(  — )) - 

L  £r[al  2 

4/ra^/3n/  ,  ,  ,n7i\,  -  e"7"* 


A2„  =  —  2^  —  (cos(a„a:)  -  cos(  —  ))- 


L  z — '  a, 

n—1 


(B.8) 

(B-9) 

(B.10) 


and  a  = 


a 

h( X  T  2p) 


,  for  an  applied  displacement  given  by  uv  = 


at  ,t  <t i 
at i  , t  >  t± 
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Figure  captions 

Fig.  1  A  rectangular  homogeneous  tissue  sample  in  an  unconfined  compression  test. 
The  sample  has  dimensions  of  t  x  ft  and  fluid  can  flow  freely  across  the  lateral 
boundaries  (p  =  0),  which  are  also  traction  free.  The  fluid  cannot  flow  across 
the  top  and  bottom  boundaries,  which  are  also  shear  stress  free  (slip  boundary 
conditions).  uy  =  0  at  the  bottom  while  a  displacement  step  (or  ramp)  function 
is  applied  at  the  top. 

Fig.  2  The  solution  for  Equation  (B.7),  at  t  =  0.3 sec,  t  =  5.4sec  and  t  =  9.9sec,  for 
L  =  10cm,  h  =  10cm,  a  =  and  t\  =  0.3sec.  For  the  blue  line,  the  poroelastic 
parameters  are  chosen  such  that  the  fluid  exchange  between  compartments  is  the 
dominant  phenomenon.  These  parameters  are  listed  in  Table  I  and  correspond  to 
the  inclusion  in  Experiment  1.  For  the  red  line,  we  increase  k  10000  times  and 
decrease  \  1000  times. 

Fig.  3  The  strain  field,  exx,  at  t  =  5.4sec,  for  both  cases  shown  in  Figure  3.  In  (a),  in 
agreement  with  the  pressure  field,  we  see  that  the  strain  is  approximately  constant 
along  almost  the  whole  sample,  indicating  that  the  percolation  is  negligible  outside 
the  very  thin  boundary  layers  located  at  the  sample  laterals.  In  (b)  we  see  that 
the  strain  varies  with  x,  indicating  that  percolation  is  important. 

Fig.  4  Experiment  1.  The  circular  inclusion  has  1cm  diameter  and  the  sample  has 
dimensions  10cm.  x  10cm.  The  fluid  cannot  flow  across  the  boundaries.  The  tissue 
is  fixed  at  the  bottom,  where  ux=uy=0.  The  lateral  surfaces  are  traction  free.  At 
the  top,  we  simulate  the  mechanical  loading  from  a  compressor  of  5 cm  of  width. 
The  displacement  of  the  compressor  is  modeled  by  a  ramp  function  such  that 
the  prescribed  uy  goes  from  0  to  1cm.  in  0.3sec,  in  the  region  corresponding  to 
x  =  2.5 cm  to  x  =  7. 5cm.  Bellow  the  compressor  we  prescribe  zero  shear  stress 
(' Tyx )• 

Fig.  5  Experiment  2.  The  circular  inclusion  has  1cm  diameter  and  the  sample  has 
dimensions  of  10cm  x  10cm.  The  fluid  cannot  flow  across  the  boundaries  and  the 
tissue  is  fixed  at  the  bottom.  The  model  is  completely  confined  at  the  top,  where 
uy  goes  from  0  to  0.03cm  in  0.3sec,  and  is  partially  confined  on  the  sides  from 
y  =  2.0 cm  to  y  =  10 cm,  while  it  is  traction  free  from  y  =  0.0cm  to  y  =  2.0 cm. 
The  goal  is  to  reproduce  a  situation  of  partial  confined  test. 

Fig.  6  (a)  The  pressure  field  (KPa)  at  t  =  0.3sec.  (b)  The  pressure  field  (KPa) 
at  t  5.4sec.  (c)  Cxx(i=5.4sec)  ^xx(i=0.3sec)  ■  (d)  eyy(t=5AseC)  ^yy(t=0.3sec)' 

(^)  (^xx(i=5. 4sec)  A  ^yy(t— 5.4sec))  xx(t=o.3sec )  T  £yy(t=o.3sec)') •  (f)  The  pressure 

field  (KPa)  at  t  10.2sec.  (g)  € xx(t—13.0sec )  ^xx(t—10.2sec)  •  (h)  ^ yy{t=13.0sec ) 

€yy(t=10.2sec)  •  (1)  xx(t—lh.Osec )  T  ^yy{t— 13.0  sec) )  (^xx(i— 10.2sec)  A  Cyy(i=10.2 sec))* 

Comparing  (a)  and  (b)  shows  that  the  inclusion  relaxes  much  faster  than 
the  surrounding  tissue  due  to  the  higher  value  of  the  microvascular  filtration 
coefficient.  In  (e)  we  see  that  in  the  first  5  seconds  the  inclusion  shrinks  as  it 
relaxes,  while  the  volume  of  the  surrounding  tissue  remains  almost  unchanged 
(dilatation  s=s  0).  In  (i)  we  see  that  from  t,  =  10.2 sec  to  t  =  15.0 sec  both  the 
inclusion  and  the  surrounding  tissue  relax  by  similar  amounts. 

Fig.  7  (a)  The  pressure  field  (KPa)  at  t  =  0.3 sec.  (b)  The  pressure  field  (KPa) 
at  t  6.3sec.  (c)  Cxx(i=6.3sec)  ^xx(t=0.3sec)  *  (d)  Cyy(f=6.3sec)  €yy(t—0.3sec)’ 

(^)  (^xx(i=6.3sec)  A  eyy(t— 5 _3sec  ))  -  (€xx(t=o.3sec)  +  tyy(t=o.3sec)) ’  (f)  The  pressure 

field  (KPa)  at  t  —  13.5sec.  (g)  ^xx{t— 31. 5sec)  ^xx(t— 13. 5sec)  •  (h)  £ yy(t=31.5sec ) 
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eyy(t— I3.5sec)*  (i)  (^xx(i=31.5sec)  T  £yy(i=: 31. 5sec))  (^xx(i=U3.5sec)  F  Cyy(t=13.5sec)') • 

In  (e)  we  see  that  the  shrinking  during  the  inclusion  relaxation  is  similar  to  the 
previous  case.  In  (i)  we  see  that  from  t  =  13.5sec  to  t  =  31.5sec  the  inclusion 
swells  while  the  surrounding  tissue  shrinks.  In  fact,  the  inclusion  attains  its 
smallest  volume  at  approximately  t  =  13.5sec. 

Fig.  8  (a)  The  pressure  held  (KPa)  at  t  =  0.3sec.  (b)  The  pressure  held 

(KPa)  at  t  5.4sec.  (c)  Cxx(i=5.4sec)  ^xx(i— 0.3sec)‘  (/I)  eyy(^— 5.4sec) 

eyyft— 0.3sec)  *  (^)  (^xx(i=5.4sec)  T  ^yy(t=5.4sec) )  (^xx(t=0.3sec)  F  £yy(t— 0.3sec))-  (f) 

The  pressure  held  (KPa)  at  t  =  10.2 sec.  (g)  exx(t=i5.0sec)  -  exx(t=io.2sec)-  (h) 

^yy(t=15.0sec)  ^yy{t—\e.2sec)‘  (i)  (bxx(i=15.0sec)  T  ^yy{t— 15.0sec) )  (^xx(t=10.2sec)  F 

tyy(t=io.2sec)) •  In  (e)  we  observe  that  the  spatio-temporal  patterns  of  strain 
isn’t  perturbed  significantly  during  the  inclusion  relaxation  by  increasing  the 
interstitial  permeability  by  100  times.  We  see  in  (i)  that  the  thin  region  around 
the  inclusion  identified  in  Figure  6(i)  is  now  much  thicker,  indicating  percolation. 

Fig.  9  (a)  The  pressure  held  (KPa)  at  t  =  0.3sec.  (b)  The  pressure  held  (KPa) 
at  t  =  5.4 sec.  (c)  The  pressure  held  (KPa)  at  t  =  10.2 sec.  Due  to  the 
small  interstitial  hydraulic  conductivity,  the  effect  of  the  permeable  boundary 
is  imperceptible  inside  the  region  of  interest,  being  confined  to  a  very  thin  layer 
located  at  the  lateral  surfaces. 

Fig.  10  A  portion  of  soft  tissue.  Its  boundary  and  volume  are  represented  by  P  and 
£2,  respectively.  Its  total  volume,  £2,  is  divided  into  three  different  compartments: 
the  interstitial  compartment,  the  hemal  vascular  compartment  and  the  lymphatic 
vascular  compartment.  The  interstitial  compartment  is  itself  a  biphasic  solid- fluid 
mixture,  where  both  fluid  and  solid  phases  can  move  independently  of  each  other. 
Tab.  1  Poroelastic  parameters  used  in  simulations. 
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